Skip to content

Core Physics

Low-level seismic physics functions: moment tensors, magnitude conversions, and the Davis (1986) point source displacement model.

Davis (1986) Point Source Model

davis_point_source 

davis_point_source(x: ndarray, y: ndarray, xcen: float, ycen: float, depth: float, Mxx: float, Myy: float, Mzz: float, Mxy: float, Myz: float, Mzx: float, nu: float = DEFAULT_POISSON_RATIO, mu: float = DEFAULT_SHEAR_MODULUS_PA) -> Tuple[np.ndarray, np.ndarray, np.ndarray]

Compute surface displacement from a point moment tensor source.

Based on Davis (1986) formulation for a point source in an elastic half-space. This is the far-field approximation valid when the observation distance is much larger than the source dimension.

Coordinate system: - x: East (positive eastward) - y: North (positive northward) - z: Up (positive upward at surface, z=0) - Depth is positive downward

PARAMETER DESCRIPTION
x

Observation coordinates in meters (can be 1D or 2D arrays)

TYPE: ndarray

y

Observation coordinates in meters (can be 1D or 2D arrays)

TYPE: ndarray

xcen

Source epicenter location in meters

TYPE: float

ycen

Source epicenter location in meters

TYPE: float

depth

Source depth in meters (positive downward)

TYPE: float

Mxx

Moment tensor components in N·m (ENU convention)

TYPE: float

Myy

Moment tensor components in N·m (ENU convention)

TYPE: float

Mzz

Moment tensor components in N·m (ENU convention)

TYPE: float

Mxy

Moment tensor components in N·m (ENU convention)

TYPE: float

Myz

Moment tensor components in N·m (ENU convention)

TYPE: float

Mzx

Moment tensor components in N·m (ENU convention)

TYPE: float

nu

Poisson's ratio (default: 0.25)

TYPE: float DEFAULT: DEFAULT_POISSON_RATIO

mu

Shear modulus in Pa (default: 30 GPa)

TYPE: float DEFAULT: DEFAULT_SHEAR_MODULUS_PA

RETURNS DESCRIPTION
Ue, Un, Uz : np.ndarray

Surface displacement components in meters: - Ue: East component (positive eastward) - Un: North component (positive northward) - Uz: Vertical component (positive upward)
Notes

The point source approximation is valid when: - Source dimension << observation distance - Source dimension << depth

Best suited for small-to-moderate earthquakes (Mw < 6.5) or far-field observations.

Examples:

>>> import numpy as np
>>> x = np.linspace(-50000, 50000, 101)
>>> y = np.linspace(-50000, 50000, 101)
>>> X, Y = np.meshgrid(x, y)
>>> # Simple vertical dip-slip source
>>> Ue, Un, Uz = davis_point_source(X, Y, 0, 0, 10000,
...     Mxx=0, Myy=0, Mzz=0, Mxy=0, Myz=1e18, Mzx=0)
Source code in src/eq_insar/core/davis.py 
 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
def davis_point_source(
    x: np.ndarray,
    y: np.ndarray,
    xcen: float,
    ycen: float,
    depth: float,
    Mxx: float,
    Myy: float,
    Mzz: float,
    Mxy: float,
    Myz: float,
    Mzx: float,
    nu: float = DEFAULT_POISSON_RATIO,
    mu: float = DEFAULT_SHEAR_MODULUS_PA,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Compute surface displacement from a point moment tensor source.

    Based on Davis (1986) formulation for a point source in an elastic
    half-space. This is the far-field approximation valid when the
    observation distance is much larger than the source dimension.

    Coordinate system:
    - x: East (positive eastward)
    - y: North (positive northward)
    - z: Up (positive upward at surface, z=0)
    - Depth is positive downward

    Parameters
    ----------
    x, y : np.ndarray
        Observation coordinates in meters (can be 1D or 2D arrays)
    xcen, ycen : float
        Source epicenter location in meters
    depth : float
        Source depth in meters (positive downward)
    Mxx, Myy, Mzz, Mxy, Myz, Mzx : float
        Moment tensor components in N·m (ENU convention)
    nu : float
        Poisson's ratio (default: 0.25)
    mu : float
        Shear modulus in Pa (default: 30 GPa)

    Returns
    -------
    Ue, Un, Uz : np.ndarray
        Surface displacement components in meters:
        - Ue: East component (positive eastward)
        - Un: North component (positive northward)
        - Uz: Vertical component (positive upward)

    Notes
    -----
    The point source approximation is valid when:
    - Source dimension << observation distance
    - Source dimension << depth

    Best suited for small-to-moderate earthquakes (Mw < 6.5) or
    far-field observations.

    Examples
    --------
    >>> import numpy as np
    >>> x = np.linspace(-50000, 50000, 101)
    >>> y = np.linspace(-50000, 50000, 101)
    >>> X, Y = np.meshgrid(x, y)
    >>> # Simple vertical dip-slip source
    >>> Ue, Un, Uz = davis_point_source(X, Y, 0, 0, 10000,
    ...     Mxx=0, Myy=0, Mzz=0, Mxy=0, Myz=1e18, Mzx=0)
    """
    # Relative coordinates from source
    x1 = x - xcen  # East offset
    x2 = y - ycen  # North offset
    d = depth      # Depth (positive down)

    # Distance from source to observation points
    R2 = x1**2 + x2**2 + d**2
    R = np.sqrt(R2)

    # Protect against singularity at/near source location
    # Use minimum distance of 1% of depth
    min_R = max(depth * 0.01, 1.0)  # At least 1 meter
    R = np.maximum(R, min_R)
    R2 = R**2

    # Powers of R
    R3 = R**3
    R5 = R**5

    # Auxiliary terms involving R + depth
    Rpd = R + d      # R + depth
    Rpd2 = Rpd**2
    Rpd3 = Rpd**3

    # Coefficients for Green's functions
    alfa = (3.0 * R + d) / (R3 * Rpd3)
    beta = (2.0 * R + d) / (R3 * Rpd2)
    eta = 1.0 / (R * Rpd2)
    psi = 1.0 / (R * Rpd)

    # Pre-computed coordinate terms
    x11a = x1**2 * alfa
    x11b = x1**2 * beta
    x22a = x2**2 * alfa
    x22b = x2**2 * beta
    anu = 1.0 - 2.0 * nu  # 1 - 2ν appears frequently

    # =========================================================================
    # Green's functions for each moment tensor component
    # =========================================================================

    # ----- Mxx (East-East) contribution -----
    G11_x = x1 / 2.0 * (-1.0 / R3 + 3.0 * x1**2 / R5 + anu * (3.0 * eta - x11a))
    G11_y = x2 / 2.0 * (-1.0 / R3 + 3.0 * x1**2 / R5 + anu * (eta - x11a))
    G11_z = -0.5 * (d / R3 - 3.0 * x1**2 * d / R5 - anu * (psi - x11b))

    # ----- Myy (North-North) contribution -----
    G22_x = x1 / 2.0 * (-1.0 / R3 + 3.0 * x2**2 / R5 + anu * (eta - x22a))
    G22_y = x2 / 2.0 * (-1.0 / R3 + 3.0 * x2**2 / R5 + anu * (3.0 * eta - x22a))
    G22_z = -0.5 * (d / R3 - 3.0 * x2**2 * d / R5 - anu * (psi - x22b))

    # ----- Mzz (Up-Up) contribution -----
    G33_x = x1 / 2.0 * (3.0 * d**2 / R5 - 2.0 * nu / R3)
    G33_y = x2 / 2.0 * (3.0 * d**2 / R5 - 2.0 * nu / R3)
    G33_z = -d / 2.0 * (-3.0 * d**2 / R5 + 2.0 * nu / R3)

    # ----- Mxy (East-North) contribution -----
    G12_x = x2 * (3.0 * x1**2 / R5 + anu * (eta - x11a))
    G12_y = x1 * (3.0 * x2**2 / R5 + anu * (eta - x22a))
    G12_z = -x1 * x2 * (-3.0 * d / R5 + anu * beta)

    # ----- Mzx (Up-East) contribution -----
    G31_x = 3.0 * x1**2 * d / R5
    G31_y = 3.0 * x1 * x2 * d / R5
    G31_z = x1 * d * (3.0 * d / R5)

    # ----- Myz (North-Up) contribution -----
    G23_x = G31_y  # = 3*x1*x2*d/R5 (symmetric)
    G23_y = 3.0 * x2**2 * d / R5
    G23_z = x2 * d * (3.0 * d / R5)

    # =========================================================================
    # Scale factor and sum contributions
    # =========================================================================

    # Point source scaling: 1/(2πμ)
    # This converts moment (N·m) to displacement (m)
    scale = 1.0 / (2.0 * np.pi * mu)

    # Sum contributions from all moment tensor components
    Ue = scale * (
        Mxx * G11_x + Myy * G22_x + Mzz * G33_x +
        Mxy * G12_x + Myz * G23_x + Mzx * G31_x
    )
    Un = scale * (
        Mxx * G11_y + Myy * G22_y + Mzz * G33_y +
        Mxy * G12_y + Myz * G23_y + Mzx * G31_y
    )
    Uz = scale * (
        Mxx * G11_z + Myy * G22_z + Mzz * G33_z +
        Mxy * G12_z + Myz * G23_z + Mzx * G31_z
    )

    return Ue, Un, Uz

Moment Tensor

double_couple_moment_tensor 

double_couple_moment_tensor(strike_deg: float, dip_deg: float, rake_deg: float, M0: float) -> Tuple[float, float, float, float, float, float]

Construct double-couple moment tensor from fault geometry.

Uses Aki & Richards (2002) convention for fault orientation: - Strike: Azimuth of fault trace, clockwise from North (0-360°) When standing on the fault, the hanging wall is on the right. - Dip: Angle from horizontal to fault plane (0-90°) Measured perpendicular to strike, toward hanging wall. - Rake: Direction of slip on fault plane (-180 to 180°) Measured from strike direction in the fault plane. - 0°: Left-lateral strike-slip - 90°: Reverse/thrust (hanging wall up) - 180° or -180°: Right-lateral strike-slip - -90°: Normal fault (hanging wall down)

PARAMETER DESCRIPTION
strike_deg

Fault strike in degrees (0-360, clockwise from North)

TYPE: float

dip_deg

Fault dip in degrees (0-90, from horizontal)

TYPE: float

rake_deg

Slip rake in degrees (-180 to 180)

TYPE: float

M0

Scalar seismic moment in N·m

TYPE: float

RETURNS DESCRIPTION
Mxx, Myy, Mzz, Mxy, Myz, Mzx : float

Moment tensor components in ENU coordinates (N·m) - Mxx: East-East component - Myy: North-North component - Mzz: Up-Up component - Mxy: East-North component - Myz: North-Up component - Mzx: Up-East component
Notes

The returned moment tensor has zero trace (Mxx + Myy + Mzz = 0), which is characteristic of a pure double-couple source.

Source code in src/eq_insar/core/moment_tensor.py 
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
def double_couple_moment_tensor(
    strike_deg: float,
    dip_deg: float,
    rake_deg: float,
    M0: float,
) -> Tuple[float, float, float, float, float, float]:
    """
    Construct double-couple moment tensor from fault geometry.

    Uses Aki & Richards (2002) convention for fault orientation:
    - Strike: Azimuth of fault trace, clockwise from North (0-360°)
             When standing on the fault, the hanging wall is on the right.
    - Dip: Angle from horizontal to fault plane (0-90°)
           Measured perpendicular to strike, toward hanging wall.
    - Rake: Direction of slip on fault plane (-180 to 180°)
           Measured from strike direction in the fault plane.
           - 0°: Left-lateral strike-slip
           - 90°: Reverse/thrust (hanging wall up)
           - 180° or -180°: Right-lateral strike-slip
           - -90°: Normal fault (hanging wall down)

    Parameters
    ----------
    strike_deg : float
        Fault strike in degrees (0-360, clockwise from North)
    dip_deg : float
        Fault dip in degrees (0-90, from horizontal)
    rake_deg : float
        Slip rake in degrees (-180 to 180)
    M0 : float
        Scalar seismic moment in N·m

    Returns
    -------
    Mxx, Myy, Mzz, Mxy, Myz, Mzx : float
        Moment tensor components in ENU coordinates (N·m)
        - Mxx: East-East component
        - Myy: North-North component
        - Mzz: Up-Up component
        - Mxy: East-North component
        - Myz: North-Up component
        - Mzx: Up-East component

    Notes
    -----
    The returned moment tensor has zero trace (Mxx + Myy + Mzz = 0),
    which is characteristic of a pure double-couple source.
    """
    # Convert to radians
    strike = np.deg2rad(strike_deg)
    dip = np.deg2rad(dip_deg)
    rake = np.deg2rad(rake_deg)

    # Fault normal vector (pointing into footwall) in ENU
    # From Aki & Richards (2002) eq. 4.88-4.89
    n_e = np.sin(dip) * np.cos(strike)
    n_n = -np.sin(dip) * np.sin(strike)
    n_u = np.cos(dip)

    # Slip vector (direction of hanging wall motion relative to footwall) in ENU
    # From Aki & Richards (2002) eq. 4.88-4.89
    s_e = np.cos(rake) * np.sin(strike) - np.sin(rake) * np.cos(dip) * np.cos(strike)
    s_n = np.cos(rake) * np.cos(strike) + np.sin(rake) * np.cos(dip) * np.sin(strike)
    s_u = np.sin(rake) * np.sin(dip)

    # Normalize vectors (should already be unit vectors, but ensure numerical stability)
    n = np.array([n_e, n_n, n_u])
    s = np.array([s_e, s_n, s_u])
    n = n / (np.linalg.norm(n) + 1e-20)
    s = s / (np.linalg.norm(s) + 1e-20)

    # Double-couple moment tensor: M = M0 * (s⊗n + n⊗s)
    M = M0 * (np.outer(s, n) + np.outer(n, s))

    # Extract components (indices: 0=E, 1=N, 2=U)
    Mxx = M[0, 0]  # Mee
    Myy = M[1, 1]  # Mnn
    Mzz = M[2, 2]  # Muu
    Mxy = M[0, 1]  # Men
    Myz = M[1, 2]  # Mnu
    Mzx = M[2, 0]  # Mue

    return Mxx, Myy, Mzz, Mxy, Myz, Mzx

mw_to_m0 

mw_to_m0(Mw: float) -> float

Convert moment magnitude to scalar seismic moment.

Uses Hanks & Kanamori (1979) relation: log10(M0) = 1.5 * Mw + 9.1 (M0 in N·m)

PARAMETER DESCRIPTION
Mw

Moment magnitude

TYPE: float

RETURNS DESCRIPTION
M0

Scalar seismic moment in N·m

TYPE: float

Examples:

>>> mw_to_m0(6.0)
1.2589254117941662e+18
>>> mw_to_m0(7.0)
3.9810717055349694e+19
Source code in src/eq_insar/core/moment_tensor.py 
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
def mw_to_m0(Mw: float) -> float:
    """
    Convert moment magnitude to scalar seismic moment.

    Uses Hanks & Kanamori (1979) relation:
        log10(M0) = 1.5 * Mw + 9.1  (M0 in N·m)

    Parameters
    ----------
    Mw : float
        Moment magnitude

    Returns
    -------
    M0 : float
        Scalar seismic moment in N·m

    Examples
    --------
    >>> mw_to_m0(6.0)
    1.2589254117941662e+18
    >>> mw_to_m0(7.0)
    3.9810717055349694e+19
    """
    return 10.0 ** (1.5 * Mw + 9.1)

m0_to_mw 

m0_to_mw(M0: float) -> float

Convert scalar seismic moment to moment magnitude.

Inverse of mw_to_m0.

PARAMETER DESCRIPTION
M0

Scalar seismic moment in N·m

TYPE: float

RETURNS DESCRIPTION
Mw

Moment magnitude

TYPE: float

Source code in src/eq_insar/core/moment_tensor.py 
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def m0_to_mw(M0: float) -> float:
    """
    Convert scalar seismic moment to moment magnitude.

    Inverse of mw_to_m0.

    Parameters
    ----------
    M0 : float
        Scalar seismic moment in N·m

    Returns
    -------
    Mw : float
        Moment magnitude
    """
    return (np.log10(M0) - 9.1) / 1.5

slip_from_moment 

slip_from_moment(M0: float, length_m: float, width_m: float, mu: float = DEFAULT_SHEAR_MODULUS_PA) -> float

Calculate average fault slip from seismic moment and fault dimensions.

Uses the relation: M0 = μ * A * D where μ is shear modulus, A is fault area, D is average slip.

PARAMETER DESCRIPTION
M0

Scalar seismic moment in N·m

TYPE: float

length_m

Fault length in meters

TYPE: float

width_m

Fault width (down-dip) in meters

TYPE: float

mu

Shear modulus in Pa (default: 30 GPa)

TYPE: float DEFAULT: DEFAULT_SHEAR_MODULUS_PA

RETURNS DESCRIPTION
slip_m

Average slip in meters

TYPE: float

Source code in src/eq_insar/core/moment_tensor.py 
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
def slip_from_moment(
    M0: float,
    length_m: float,
    width_m: float,
    mu: float = DEFAULT_SHEAR_MODULUS_PA,
) -> float:
    """
    Calculate average fault slip from seismic moment and fault dimensions.

    Uses the relation: M0 = μ * A * D
    where μ is shear modulus, A is fault area, D is average slip.

    Parameters
    ----------
    M0 : float
        Scalar seismic moment in N·m
    length_m : float
        Fault length in meters
    width_m : float
        Fault width (down-dip) in meters
    mu : float
        Shear modulus in Pa (default: 30 GPa)

    Returns
    -------
    slip_m : float
        Average slip in meters
    """
    area = length_m * width_m
    return M0 / (mu * area)