Seismology

Seismology: earthquakes on loaded faults.

FaultElement(label, visual=None, *, force=0.0)

Bases: Element

Block of the Olami-Feder-Christensen earthquake model, the force on it.

The cascade of earthquakes loads every block uniformly until one reaches the threshold force, and a slipping block passes a fraction of its force to each of its four neighbours and drops to zero.

Parameters:

Name	Type	Description	Default
label	str	The label of the block or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
force	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial force.	0.0

Examples:

>>> from skcomplex.models.physics.seismology import FaultElement
>>> FaultElement('block').force
0.0

Source code in src/skcomplex/models/physics/seismology/_faults.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    force: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.force = force

SpringBlockElement(label, visual=None, *, x=0.0, v=0.0, m=1.0, k_c=2.5, k_p=1.0, velocity=0.01, F0=1.0, v_c=0.5)

Bases: Element

Block of the Burridge-Knopoff fault, pulled by springs and a moving plate against a friction that weakens.

The force on a block is the pull k_c of the springs to its neighbours plus the pull k_p of the plate moving at velocity. A block at rest sticks until the force exceeds the static friction F0, and a sliding block feels the friction F0 / (1 + v / v_c). The stick condition reads the whole force, the coupling included, so the block keeps its own rate laws.

Parameters:

Name	Type	Description	Default
label	str	The label of the block or the field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial displacement.	0.0
v	Annotated[float, StateSpec()]	The initial velocity.	0.0
m	Annotated[float, ParameterSpec()]	The mass.	1.0
k_c	Annotated[float, ParameterSpec()]	The stiffness of the springs between neighbouring blocks.	2.5
k_p	Annotated[float, ParameterSpec()]	The stiffness of the springs to the plate.	1.0
velocity	Annotated[float, ParameterSpec()]	The speed of the plate.	0.01
F0	Annotated[float, ParameterSpec()]	The static friction.	1.0
v_c	Annotated[float, ParameterSpec()]	The characteristic speed of the weakening of the friction.	0.5

Examples:

>>> from skcomplex.models.physics.seismology import SpringBlockElement
>>> SpringBlockElement('block').F0
1.0

Source code in src/skcomplex/models/physics/seismology/_spring_blocks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    v: Annotated[float, StateSpec()] = 0.0,
    m: Annotated[float, ParameterSpec()] = 1.0,
    k_c: Annotated[float, ParameterSpec()] = 2.5,
    k_p: Annotated[float, ParameterSpec()] = 1.0,
    velocity: Annotated[float, ParameterSpec()] = 0.01,
    F0: Annotated[float, ParameterSpec()] = 1.0,
    v_c: Annotated[float, ParameterSpec()] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.v = v
    self.m = m
    self.k_c = k_c
    self.k_p = k_p
    self.velocity = velocity
    self.F0 = F0
    self.v_c = v_c

burridge_knopoff(n_blocks=20, *, k_c=2.5, k_p=1.0, velocity=0.01, F0=1.0, v_c=0.5, random_state=0)

Return the Burridge-Knopoff fault, a chain of blocks pulled by a moving plate that sticks and slips.

Every block is coupled to its neighbours by springs of stiffness k_c and to a plate moving at velocity by a leaf spring of stiffness k_p. A block at rest stays stuck until its spring force exceeds the static friction F0, and a sliding block feels the friction F0 / (1 + v / v_c), which weakens with its speed, so slips spread as earthquakes of many sizes. The ends of the chain are free, and the blocks start at small random offsets.

The drift switches whenever a block sticks or starts to slide, and these switches of single blocks give no events, so the adaptive default solver stalls on them: LSODA gives way to RK45 after its stall guard fires, and 2000 time units take about 24 s. The registered solver is RK4Solver(dt=0.05), well inside the stability limit set by the highest frequency sqrt(k_p + 4 k_c), which runs the same span in about 3.5 s with the same statistics of slips.

Parameters:

Name	Type	Description	Default
n_blocks	int	The number of blocks.	20
k_c	float	The stiffness of the springs between neighbouring blocks.	2.5
k_p	float	The stiffness of the leaf springs to the plate.	1.0
velocity	float	The velocity of the plate.	0.01
F0	float	The static friction force.	1.0
v_c	float	The speed at which the friction halves.	0.5
random_state	int | Generator | None	The seed or generator of the initial offsets.	0

Returns:

Name	Type	Description
system	System	The model, with the positions fault.x and velocities fault.v of the blocks.

Examples:

>>> from skcomplex.models.physics.seismology import burridge_knopoff
>>> burridge_knopoff(5).elements[0].check().prepare().size_
5

Source code in src/skcomplex/models/physics/seismology/_spring_blocks.py

@register_model(
    field='physics/seismology',
    tags=('field', 'stick-slip', 'earthquakes', 'friction'),
    references=('Burridge and Knopoff 1967', 'Carlson and Langer 1989', 'Charbonneau 2017'),
    defaults={'until': 200.0, 'dt': 1.0, 'solver': RK4Solver(dt=0.05)},
)
def burridge_knopoff(
    n_blocks: int = 20,
    *,
    k_c: float = 2.5,
    k_p: float = 1.0,
    velocity: float = 0.01,
    F0: float = 1.0,
    v_c: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Burridge-Knopoff fault, a chain of blocks pulled by a moving plate that sticks and slips.

    Every block is coupled to its neighbours by springs of stiffness `k_c` and to a plate moving at `velocity` by a
    leaf spring of stiffness `k_p`. A block at rest stays stuck until its spring force exceeds the static friction `F0`,
    and a sliding block feels the friction `F0 / (1 + v / v_c)`, which weakens with its speed, so slips spread as
    earthquakes of many sizes. The ends of the chain are free, and the blocks start at small random offsets.

    The drift switches whenever a block sticks or starts to slide, and these switches of single blocks give no events,
    so the adaptive default solver stalls on them: LSODA gives way to RK45 after its stall guard fires, and 2000 time
    units take about 24 s. The registered solver is `RK4Solver(dt=0.05)`, well inside the stability limit set by the
    highest frequency `sqrt(k_p + 4 k_c)`, which runs the same span in about 3.5 s with the same statistics of slips.

    Args:
        n_blocks:
            The number of blocks.

        k_c:
            The stiffness of the springs between neighbouring blocks.

        k_p:
            The stiffness of the leaf springs to the plate.

        velocity:
            The velocity of the plate.

        F0:
            The static friction force.

        v_c:
            The speed at which the friction halves.

        random_state:
            The seed or generator of the initial offsets.

    Returns:
        system:
            The model, with the positions `fault.x` and velocities `fault.v` of the blocks.

    Examples:
        >>> from skcomplex.models.physics.seismology import burridge_knopoff
        >>> burridge_knopoff(5).elements[0].check().prepare().size_
        5
    """
    offsets = np.random.default_rng(random_state).uniform(-0.5, 0.5, n_blocks)
    fault = FieldElement(
        'fault',
        element=SpringBlockElement,
        space=LatticeSpace(n_blocks, boundary='reflecting'),
        boundary={'x': 'zero-flux', 'v': 'zero-flux'},
        values={'x': offsets, 'k_c': k_c, 'k_p': k_p, 'velocity': velocity, 'F0': F0, 'v_c': v_c},
    )
    return System(elements=[fault])

earthquakes(shape=(50, 50), *, alpha=0.2, random_state=0)

Return the Olami-Feder-Christensen earthquake model with open edges, started from random forces.

The fault is a cascade that resets a slipping block to zero and passes the fraction alpha of its force to each of its four neighbours. With alpha below 1/4 the model is dissipative and its earthquake sizes still follow a power law. Its earthquakes are recorded as earthquake events with their sizes, the number of slips, their durations, and their areas.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(50, 50)
alpha	float	The fraction of force passed to each neighbour.	0.2
random_state	int | Generator | None	The seed or generator of the initial forces.	0

Returns:

Name	Type	Description
system	System	The model, with the forces fault.force.

Examples:

>>> from skcomplex.models.physics.seismology import earthquakes
>>> earthquakes((10, 10), alpha=0.15).interactions[0].share
0.15

Source code in src/skcomplex/models/physics/seismology/_faults.py

@register_model(
    field='physics/seismology',
    tags=('lattice', 'earthquakes', 'power law'),
    references=('Olami, Feder, and Christensen 1992', 'Jensen 2023'),
    defaults={'until': 5000, 'dt': 100},
)
def earthquakes(
    shape: Sequence[int] = (50, 50),
    *,
    alpha: float = 0.2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Olami-Feder-Christensen earthquake model with open edges, started from random forces.

    The fault is a cascade that resets a slipping block to zero and passes the fraction `alpha` of its force to each of
    its four neighbours. With `alpha` below 1/4 the model is dissipative and its earthquake sizes still follow a power
    law. Its earthquakes are recorded as `earthquake` events with their sizes, the number of slips, their durations,
    and their areas.

    Args:
        shape:
            The number of rows and columns.

        alpha:
            The fraction of force passed to each neighbour.

        random_state:
            The seed or generator of the initial forces.

    Returns:
        system:
            The model, with the forces `fault.force`.

    Examples:
        >>> from skcomplex.models.physics.seismology import earthquakes
        >>> earthquakes((10, 10), alpha=0.15).interactions[0].share
        0.15
    """
    forces = np.random.default_rng(random_state).random(tuple(shape))
    space = LatticeSpace(tuple(shape), boundary='fixed')
    return System(
        elements=[LatticeElement('fault', element=FaultElement, space=space, values={'force': forces})],
        interactions=[
            CascadeInteraction(
                'earthquake',
                target='fault',
                state='force',
                threshold=THRESHOLD,
                share=alpha,
                mode='reset',
                drive=_load,
                neighborhood=VonNeumannNeighborhood(center=True),
            ),
        ],
    )