Four tools for seismic data

Originally published at amansinghnp.com.

Over 2025 and 2026 I built four small open-source tools that sit at different points of a seismic data pipeline. Each solves a different problem, so start with the one you need:

Each tool's section opens with a short entry: the clearest result, the main limitation and how it was tested. The physics, numerics and full test results follow. Every chart below was produced by running the published code (or, where noted, a line-for-line NumPy port of it) on the data listed in section 6. Some of the results are flattering and some are not. I've kept both. The earthquake-severity model that led to seismic-linter is on the Projects page as South Asia Earthquake Severity.

< 0.05% friction parameter error after inversion

−20.7 dB ground roll removed by the f-k mute

≈ 2 dB less signal needed by x² STA/LTA at a 1% false-alarm rate

1. Why build a toolchain

Machine learning in seismology has moved fast. Deep networks now pick phases and detect events at scale (Mousavi & Beroza, 2022), and physics-informed neural networks (PINNs) put governing equations into the loss function (Raissi et al., 2019; Karniadakis et al., 2021). Two problems keep coming up anyway. The first is that the physics is often bolted on after the fact, instead of living inside the computational graph where gradients can reach it. The second is that evaluation is fragile: time-series leakage has inflated reported results across many fields (Kapoor & Narayanan, 2023), and earthquake forecasting has its own history of complex models that turn out no better than simple baselines (Mignan & Broccardo, 2019, 2020).

The four tools line up against those two problems. torch-friction and Basalt handle the physics and signal-processing side. fast-trigger handles detection. seismic-linter guards the evaluation.

2. torch-friction: differentiable rate-and-state friction

torch-friction in brief

• Result. In a teacher–student test at a lab-scale D_c of 10 µm (Setup B, section 2.5), gradient descent through the solver recovered a, b and D_c to within 0.05% from starting guesses that were 33–100% wrong, and the loss levelled off at the noise variance.
• Main limitation. Recovery depends on the experiment, not the optimiser. With only 8% of D_c of slip after a single step (Setup A), b and D_c can't be separated and b moved further from the truth. The layer also solves only the zero-dimensional friction ODE, not the elastodynamic problem around a fault.
• How it was tested. Synthetic velocity-step data from a teacher layer with hidden parameters and Gaussian noise (σ = 10⁻⁴), fitted by a student layer with Adam; RK4 checked against the exact aging-law relaxation at seven step sizes (error falls 19.7× then 16.8× per halving of Δt). Run with torch-friction from PyPI in float64 on CPU.

2.1 The friction law

Laboratory rock friction isn't a constant. It depends on how fast the surfaces slide and on how long the contacts have had to age. Dieterich (1979) and Ruina (1983) captured this with a friction coefficient that depends on slip velocity V and a state variable θ, which you can read as the average lifetime of the asperity contacts:

μ(V, θ) = μ₀ + a ln(V / V₀) + b ln(V₀θ / D_c)

The a term is the direct effect: friction jumps immediately when velocity changes. The b term is the evolution effect: as the contact population turns over across a characteristic slip distance D_c, friction relaxes to a new steady state. At steady state, θ = D_c/V and μ_ss = μ₀ + (a − b) ln(V/V₀). When a − b < 0 the fault is velocity-weakening: sliding faster makes it weaker, which is the precondition for earthquake nucleation (Marone, 1998; Scholz, 1998).

Whether that weakening actually produces an earthquake depends on the surroundings. In the classic spring-slider model, a block dragged through a spring of stiffness k slides stably if k is above a critical stiffness k_c = σ(b − a)/D_c, and stick-slips below it (Ruina, 1983; Rice & Ruina, 1983).

Figure 1. The spring-slider model. The load point moves at a constant rate, stress builds in the spring while the block sticks, then drops when the block slips. With velocity-weakening friction and k < k_c, this cycle repeats indefinitely: the laboratory analogue of the seismic cycle.

2.2 Two ways the state can evolve

The friction equation says nothing about how θ changes over time. There are two standard choices, and torch-friction implements both as law='aging' and law='slip':

The standard lab test is a velocity step. Figure 2 runs one through FrictionLayer.forward_trajectory with a = 0.015, b = 0.020 and D_c = 10 µm, stepping from 1 to 10 µm/s and back.

Figure 2. Velocity step up at t = 2 s and back down at t = 8 s, computed with torch-friction (Δt = 2 ms, RK4). The jump of a·ln10 ≈ 0.035 is the direct effect; the settling to (a − b)·ln10 ≈ −0.012 below the starting level is the evolution effect. The two laws agree on the peak and the steady state but differ on the path between them.

The two laws only differ in the transient, but the difference isn't small: up to Δμ = 0.012 about one second after the step up, which is a third of the whole direct effect. The asymmetry is easy to see from the equations. After a tenfold step up, Vθ/D_c = 10, so the aging law evolves at a rate of 1 − 10 = −9 while the slip law evolves at −10 ln 10 ≈ −23. The slip law weakens faster. After the step down (Vθ/D_c = 0.1) it's reversed: 0.9 for aging against 0.23 for slip, so the aging law heals faster. Six seconds after the step down, the aging curve has recovered to μ = 0.586 while the slip curve is still at 0.575. Lab experiments with large velocity steps tend to favour the slip law (Bhattacharya et al., 2015), which is why the library keeps both.

2.3 Integrating the state equation

Each call to the layer advances θ by one classical fourth-order Runge-Kutta (RK4) step, then evaluates the friction equation at the new state. The core is short:

torch_friction/layers.py

def rk4_step(self, t, theta, v, dt, Dc):
    k1 = self.physics(t, theta, v, Dc)
    k2 = self.physics(t + dt/2, theta + k1*dt/2, v, Dc)
    k3 = self.physics(t + dt/2, theta + k2*dt/2, v, Dc)
    k4 = self.physics(t + dt, theta + k3*dt, v, Dc)
    return theta + (dt/6.0)*(k1 + 2*k2 + 2*k3 + k4)

The state equation is stiff when D_c is small relative to V·Δt. Linearising the aging law about steady state gives an eigenvalue of −V/D_c, and RK4's stability interval on the negative real axis ends near −2.79 (Hairer & Wanner, 1996), so explicit integration needs Δt < 2.79 D_c/V. The layer warns once Δt goes past 0.5 D_c/V_max, which leaves a margin of about 5.6× below the stability limit. That margin is sensible, because the accuracy limit arrives well before the stability limit does.

To check the accuracy, Figure 3 relaxes the aging law from θ₀ = 10 s at constant V, where the exact solution is θ(t) = D_c/V + (θ₀ − D_c/V)e^−Vt/D_c, and compares it with the library's RK4 at seven step sizes.

Figure 3. Error in θ after 2 s of relaxation against step size. Halving Δt cuts the error by 19.7× and then 16.8×, close to the 2⁴ = 16 of a fourth-order method. Below Δt ≈ 0.06 s the error stops falling at about 3 × 10⁻⁷ s.

The floor isn't RK4's fault. The layer guards its logarithms and divisions with a 10⁻¹² epsilon and a smooth soft_clamp, and with D_c = 10⁻⁵ m those guards shift the steady state by about one part in 10⁷. That's a deliberate trade: gradients that never hit NaN in exchange for a relative bias far below lab measurement noise. If you're working at very small D_c, scale the epsilon with it.

2.4 Backpropagating through the solver

Because every operation in the RK4 step is a PyTorch tensor op, autograd records the whole integration. forward_trajectory unrolls the step T times into preallocated output tensors, and loss.backward() replays that graph in reverse to get ∂L/∂a, ∂L/∂b, ∂L/∂D_c and ∂L/∂μ₀ (Paszke et al., 2019). This is "discretise-then-optimise": the gradients are exact for the discrete solver, at the cost of memory that grows linearly with T. The alternative is the adjoint method used by neural ODEs, which solves a second ODE backwards in time and uses constant memory, but its gradients are only as accurate as that backward solve (Chen et al., 2018). For sequences of a few thousand steps, exact gradients are the better deal.

Figure 4. One cell of the unrolled graph. Black arrows are the forward pass; green arrows are gradients flowing back through all four RK4 stages to the parameters. Parameters are stored in unconstrained form and passed through softplus, so a, b and D_c stay positive whatever the optimiser does.

Three other details matter for training. Parameters default to float64, because the evolution term takes logs of quantities that span many orders of magnitude. Negative velocities raise an error instead of being silently clamped. And freeze_parameters(['a', 'Dc']) lets you fix quantities you've measured independently, which matters a great deal for the next section.

2.5 The inverse problem, and when it's solvable

The test is to generate noisy friction data from a "teacher" layer with hidden parameters, then train a "student" layer, started from wrong values, to recover them. I ran two versions.

Setup A is the demo in the repository's examples/ folder: one step from 1 to 10 µm/s, D_c = 10 mm, Δt = 0.1 s, 1,000 samples, Adam at a learning rate of 10⁻³ for 150 epochs. Setup B uses a lab-scale D_c = 10 µm, a sequence of three steps (up ×10, down ×10, up ×3) over 12 s at Δt = 2 ms, and Adam at 10⁻² for 600 epochs. Both add Gaussian noise with σ = 10⁻⁴.

Figure 5. Setup B: relative error of each recovered parameter against epoch, starting from guesses that were 33–100% wrong. All three are within 1% by epoch 166 and within 0.05% at the end. The loss levels off at 9.9 × 10⁻⁹, which is the noise variance (10⁻⁸): the student is fitting everything except the noise.

Setup A doesn't recover the parameters, and b actually moves further from the truth. Part of the reason is optimisation (the loss is still 1,000× the noise floor after 150 epochs), but the main reason is identifiability. After the step, the block slides at 10 µm/s for 80 s, which is 0.8 mm, only 8% of D_c. The evolution effect has barely started, so the data mostly constrain the direct effect plus a linearised evolution term in which b and D_c appear together as roughly b/D_c. Many (b, D_c) pairs fit equally well. Setup B slides many times D_c after each step and steps in both directions, so the two parameters leave separate fingerprints.

Design rule

For rate-and-state inversion, design the experiment before the optimiser. Record at least several D_c of slip after each step, include steps up and down, and freeze any parameter you've measured some other way. A differentiable solver finds the minimum; it can't create information the data don't contain.

2.6 Where this sits in the literature

torch-friction is deliberately narrow. It solves the zero-dimensional friction ODE, not the elastodynamic problem around a fault. Rucker and Erickson (2024) go further and train PINNs for the wave equation with a rate-and-state fault boundary in 1D and 2D, recovering depth-varying friction parameters. Borate et al. (2024) show a PINN can recover rate-and-state parameters from acoustic monitoring of lab stick-slip experiments. A differentiable friction layer is a building block for that kind of model: a component whose physics is exact, placed inside a larger network that learns everything else.

3. Basalt: f-k filtering at the metal

Basalt in brief

• Result. On a synthetic shot gather, the sigmoid mask with v_cut = 1,500 m/s kept 99.97% of the reflection energy and 0.85% of the ground-roll energy, a 20.7 dB reduction (section 3.1).
• Main limitation. There are no throughput numbers yet: every results cell in PERFORMANCE.md is marked TBD, and the C++ AVX2 kernels weren't run for this post. The kernel also reads only mute_vel_max, so the reject cone is always symmetric in ±k.
• How it was tested. A NumPy port of Basalt's mask applied to a synthetic gather with separately known reflection and ground-roll components, with energy measured on each component. In the repository's CI, a golden checksum for single-thread AVX2 output and a bit-for-bit match between multi-threaded batch and single-threaded output.

3.1 What an f-k filter does

In an active-source shot gather, the reflections geologists care about arrive with high apparent velocities across the receiver spread. Ground roll, the surface waves travelling along the ground, is slow (typically a few hundred m/s), low-frequency and much stronger. A 2D Fourier transform maps each linear event with apparent velocity v to a line f = v·k in the frequency-wavenumber plane, so the two separate by slope. Muting the slow wedge and transforming back is the "pie-slice" filter (Embree et al., 1963; Yilmaz, 2001).

Basalt's mute is a smooth sigmoid rather than a hard cut, because a hard edge in the f-k plane rings in the time domain:

w(f, k) = sigmoid( 10 · ( |f / (k + ε)| − v_cut ) / taper )

Figure 6. (a) A synthetic gather: 128 traces at 5 m spacing, 4 ms sampling, four reflectors (2,000–3,600 m/s, 30 Hz Ricker) and two ground-roll modes (320 m/s at 8 Hz, 450 m/s at 12 Hz). (b) Its f-k spectrum; the shaded region, |f/k| < 1,500 m/s, is muted. (c) The result, using a NumPy port of Basalt's mask with v_cut = 1,500 m/s and taper = 100 m/s.

Measured on each component separately (the filter is linear), the mask keeps 99.97% of the reflection energy and 0.85% of the ground-roll energy, a 20.7 dB reduction. The panel also shows the two classic f-k artefacts. The small residue at the top centre is the ground-roll apex at zero offset, where the event is locally flat, so its energy sits near k = 0 inside the pass region. The faint smears at the edges near 1.0 s are wrap-around: the FFT treats the gather as periodic in offset, so energy that gets cut at one edge leaks back in at the other. Production workflows handle these with spatial tapering and padding before the transform.

3.2 Pipeline and memory

Figure 7. Basalt's path from disk to filtered data. The lower strip shows the mask kernel sweeping one frequency row eight wavenumbers at a time, one AVX2 register per iteration.

Most of the engineering is about not moving data. SEG-Y files are memory-mapped rather than read into buffers, and the parser decodes the EBCDIC textual header and legacy IBM-360 floating-point samples in place (SEG Technical Standards Committee, 2017). All working memory comes from a MemoryArena, a bump allocator with 64-byte alignment, so one gather's scratch space is a single allocation and reset() frees it in O(1). Complex spectra are stored as a structure of arrays (separate real and imag planes), so that eight consecutive real parts load into one 256-bit register without shuffling.

A 2D FFT is a set of 1D FFTs along rows, a transpose, then 1D FFTs along what used to be columns (Cooley & Tukey, 1965). The transpose is where naive code falls apart: striding down columns of a large matrix misses the cache on nearly every access. Basalt tiles the transpose so each block fits in cache, the approach cache-oblivious algorithms formalise (Frigo et al., 1999), and picks the 2D FFT schedule automatically by problem size, in the same spirit as FFTW's planner (Frigo & Johnson, 2005).

3.3 The SIMD mask kernel

The filter itself is embarrassingly parallel: every f-k bin gets multiplied by its own weight. The AVX2 path tiles the spectrum (32×64 bins for small gathers, 64×128 by default, 32×96 above 2,048), prefetches 8 to 24 elements ahead, and handles eight bins per instruction:

src/filter/fk_filter.cpp (AVX2 path, abridged)

const __m256 vk     = _mm256_load_ps(k_vals.data());
const __m256 vel    = _mm256_div_ps(vf, _mm256_add_ps(vk, v_eps));
__m256 x            = _mm256_sub_ps(absolute_ps(vel), v_cut);
x                   = _mm256_mul_ps(x, v_steep);
const __m256 weight = fast_sigmoid_avx(x);
real = _mm256_mul_ps(real, weight);
imag = _mm256_mul_ps(imag, weight);

There are no branches in the inner loop: the "mute or keep" decision is a multiplication by a weight between 0 and 1, so no lane waits on a mispredicted branch. At runtime, the library uses the AVX-512 path if the CPU reports avx512f, vl, bw and dq, and otherwise falls back to AVX2. One detail worth knowing: the kernel reads only mute_vel_max, so the reject cone is always symmetric in ±k. mute_vel_min is currently unused.

3.4 Performance: what's measured and what isn't

I'm not quoting throughput numbers here, because the repository hasn't published any yet. PERFORMANCE.md lays out the method: pinned CPU governors, the basalt_bench sweep across size, threads, SIMD mode, NUMA mode and memory policy, and PMU counters on self-hosted Linux. Every results cell is still marked TBD, because shared CI runners are fine for correctness testing but no good for performance claims. What is in CI is correctness: a golden checksum for single-thread AVX2 output, and a check that multi-threaded batch output matches single-threaded output bit for bit.

4. fast-trigger: STA/LTA in C++23

fast-trigger in brief

• Result. The detector triggers on both P and S in a synthetic seismogram. The sharper finding is about the characteristic function: at a fixed 1% false-trigger rate, ObsPy's x² version reaches 50% detection at about 1.8σ against about 2.25σ for the |x| that fast-trigger uses, roughly 2 dB less signal, and picks earlier (a median 0.26 s after onset at 3σ, against 0.39 s). See section 4.3.
• Main limitation. The characteristic function is fixed at |x|, which the test above shows is the weaker detector for this signal and noise model. It hasn't yet been benchmarked on real waveforms, and the package isn't on PyPI yet.
• How it was tested. An exact NumPy port of the recursion, compared with ObsPy 1.5.1's recursive_sta_lta on synthetic 100 Hz seismograms with known P and S onsets. Thresholds were calibrated on 300 noise-only traces, with 150 trials per P-wave amplitude and fixed seeds. The repository adds Hypothesis property-based tests and an AddressSanitizer and UBSan build.

4.1 Classic and recursive STA/LTA

The short-term-average over long-term-average ratio is still the most widely used automatic event detector (Allen, 1978; Trnkoczy, 2012). A short window tracks the signal's recent energy, a long window tracks the background, and a trigger fires when their ratio crosses a threshold. The recursive form replaces the sliding windows with exponential moving averages:

STA_i = c_s|x_i| + (1 − c_s)STA_i−1,    c_s = 1/N_sta

Withers et al. (1998) compare these and other trigger algorithms in detail. fast-trigger implements the recursive form in C++23, starts both averages at the first sample's magnitude, and exposes compute_sta_lta for a single trace and compute_sta_lta_batch for a 2D array of traces, the latter parallelised over traces with OpenMP. Inputs are validated at the boundary: empty arrays, zero windows, STA ≥ LTA and windows longer than the trace all raise ValueError in Python.

4.2 A worked example

Figure 8 applies it to a synthetic 100 Hz seismogram: Gaussian noise plus a 0.2 Hz microseism, a 6 Hz P arrival at 30 s and a stronger 3 Hz S arrival at 40 s. STA is 0.5 s and LTA is 10 s. For comparison, it also shows ObsPy's recursive_sta_lta (Beyreuther et al., 2010), which uses the same recursion on x² instead of |x|.

Figure 8. Top: the seismogram (min/max envelope). Bottom: STA/LTA from both characteristic functions, with each one's threshold set for a 1% false-trigger rate on noise alone (dashed). Both trigger on P and S; the |x| ratio triggers at 30.23 s and the x² ratio at 30.16 s. The green line sweeps the trace the way a real-time detector would.

4.3 The characteristic function matters more than the language

Squaring the signal doubles its dynamic range in decibels. Where an arrival lifts the |x| ratio to about 3.4, it lifts the x² ratio to about 7.7. The noise-only ratio gets peakier too, though, so the only fair comparison is at an equal false-alarm rate. I calibrated both thresholds on 300 noise-only traces (1.50 for |x|, 1.92 for x²), then measured how often each detects a P wave of known amplitude over 150 trials per level.

Figure 9. Detection probability against P-wave amplitude at a fixed 1% false-trigger rate. The x² function reaches 50% detection at about 1.8σ; |x| needs about 2.25σ, roughly 2 dB more signal. The x² function also picks earlier: a median 0.26 s after onset at 3σ, against 0.39 s.

For this signal and noise model, the energy function is simply the better detector, and no amount of C++ speed makes up for missing an event. Making the characteristic function a parameter (|x|, x², or Allen's (1978) function, which adds a weighted derivative term) is the most useful next change to fast-trigger. It's also a reminder to benchmark detectors at a fixed false-alarm rate, not at a fixed threshold.

4.4 Engineering notes

The C++ core takes std::span inputs and uses an mdspan-style view for batches, with __restrict pointers so the compiler can vectorise. Release builds on Linux and macOS use -ffast-math by default, which is fine for a detector but not for bit-reproducible research, so FAST_TRIGGER_STRICT_FLOAT=ON restores strict IEEE behaviour. There are Hypothesis property-based tests, and an AddressSanitizer and UBSan build. One practical caveat: because the batch function already uses OpenMP internally, calling it from many Python threads at once oversubscribes the CPU. Use one batch call at a time, or set OMP_NUM_THREADS.

5. seismic-linter: catching leakage before training

seismic-linter in brief

• Result. Version 0.2.0 flagged all three leakage mistakes in an eight-line script (four violations across rules T001, T002 and T003) and exited with code 1, which is what lets it block a merge from a pre-commit hook or a GitHub Action (section 5.2).
• Main limitation. The rules are heuristics over syntax, not data flow. T002 judges by variable names, so a leaky fit on train_df gets through; T001 doesn't yet treat ewm or custom wrappers as safe; T003 only accepts a literal shuffle=False.
• How it was tested. Shown here on one deliberately leaky script, with the real output reproduced in section 5.2. Two runtime checks, @verify_monotonicity and validate_split_integrity, cover what static analysis can't see.

5.1 Three ways the future leaks into the past

Kaufman et al. (2012) define leakage as information in the training data that wouldn't be available at prediction time. In time series it usually comes from ordinary-looking code. Normalising with a global mean includes test-period statistics in every training row. Fitting a scaler before splitting does the same thing. And a shuffled train/test split lets the model train on 2023 and be tested on 2019. Standard k-fold cross-validation carries the same risk for time series unless it's done with care (Bergmeir & Benítez, 2012). Across 17 fields, Kapoor and Narayanan (2023) found leakage affecting 294 published papers.

Figure 10. Top: a shuffled split scatters test years (red) through the record, so the model trains on years that come after the ones it's tested on. Bottom: the temporal split used in the South Asia earthquake severity project, where every validation and test event is later than every training event.

5.2 How the analysis works

Figure 11. seismic-linter's pipeline. Notebooks are flattened to code cells, each file is parsed once and cached by content hash, and files are analysed in parallel processes.

The linter walks Python's abstract syntax tree and reasons about context. A call to .mean() is fine if it follows groupby, rolling, expanding or resample in the same chain, because those carry temporal boundaries; on a bare DataFrame, it gets flagged.

Here's the real output of version 0.2.0 on an eight-line script that makes all three mistakes:

$ seismic-linter lintdemo/
File: lintdemo/train.py
  ❌ [Line 6] T001: Global mean() may cause temporal leakage.
  ❌ [Line 6] T001: Global std() may cause temporal leakage.
  ℹ️ [Line 7] T002: Model fit() called on potential non-training data 'df'.
  ❌ [Line 8] T003: train_test_split used with shuffle=True on time series

Found 4 potential violation(s).        # exit code 1

A non-zero exit code is what turns this into a gate: it can run as a pre-commit hook or a GitHub Action, and fail-on / ignore in pyproject.toml control which rules block a merge. Deliberate exceptions get an inline # seismic-linter: ignore T001.

5.3 Runtime guards and limits

Static analysis can't see data, so the package also ships two runtime checks: a @verify_monotonicity(time_column="time") decorator that rejects unsorted catalogues, and validate_split_integrity(train, test, time_column="time"), which raises TemporalCausalityError if any test timestamp is earlier than a training timestamp.

The static rules are heuristics, and they have limits. T002 decides from variable names, so a leaky fit on a variable called train_df gets through. T001 doesn't yet recognise ewm or custom wrappers as safe, so those produce false positives. T003 only accepts a literal shuffle=False, so a variable holding False still gets reported. A proper fix would be data-flow analysis that follows which rows reach each fit call. For a linter whose job is to be run on every commit, false positives you can suppress are a better failure mode than silent misses.

6. Datasets

I've used synthetics wherever the point was to test an algorithm, because only synthetic data come with exact ground truth. The two "not used yet" rows are where each tool should go next.

7. Reproducibility notes

• Friction figures use torch-friction from PyPI (PyTorch 2.14, float64, CPU). The Setup A and B scripts follow examples/inverse_problem_demo.py, with the changes described in section 2.5.
• STA/LTA figures use an exact NumPy port of fast-trigger's recursion (same coefficients, same first-sample initialisation), because the package isn't on PyPI yet, and ObsPy 1.5.1 for the x² comparison. Seeds are fixed.
• The f-k figure uses a NumPy port of Basalt's sigmoid mask. The C++ kernels need x86-64 with AVX2, so they weren't run for this post.

You can backpropagate through physics. A friction law can be a layer in a neural network.

References

1. Allen, R. V. (1978). Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68(5), 1521–1532. https://doi.org/10.1785/BSSA0680051521
2. Bergmeir, C., & Benítez, J. M. (2012). On the use of cross-validation for time series predictor evaluation. Information Sciences, 191, 192–213. https://doi.org/10.1016/j.ins.2011.12.028
3. Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y., & Wassermann, J. (2010). ObsPy: A Python toolbox for seismology. Seismological Research Letters, 81(3), 530–533. https://doi.org/10.1785/gssrl.81.3.530
4. Bhattacharya, P., Rubin, A. M., Bayart, E., Savage, H. M., & Marone, C. (2015). Critical evaluation of state evolution laws in rate and state friction: Fitting large velocity steps in simulated fault gouge with time-dependent healing. Journal of Geophysical Research: Solid Earth, 120(9), 6365–6385. https://doi.org/10.1002/2015JB012437
5. Borate, P., Rivière, J., Marty, S., Marone, C., Kifer, D., & Shokouhi, P. (2024). Physics informed neural network can retrieve rate and state friction parameters from acoustic monitoring of laboratory stick-slip experiments. Scientific Reports, 14, Article 24624. https://doi.org/10.1038/s41598-024-75826-y
6. Chen, R. T. Q., Rubanova, Y., Bettencourt, J., & Duvenaud, D. (2018). Neural ordinary differential equations. In Advances in Neural Information Processing Systems 31. https://arxiv.org/abs/1806.07366
7. Cooley, J. W., & Tukey, J. W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90), 297–301. https://doi.org/10.1090/S0025-5718-1965-0178586-1
8. Dieterich, J. H. (1979). Modeling of rock friction: 1. Experimental results and constitutive equations. Journal of Geophysical Research: Solid Earth, 84(B5), 2161–2168. https://doi.org/10.1029/JB084iB05p02161
9. Embree, P., Burg, J. P., & Backus, M. M. (1963). Wide-band velocity filtering: The pie-slice process. Geophysics, 28(6), 948–974. https://doi.org/10.1190/1.1439310
10. Frigo, M., & Johnson, S. G. (2005). The design and implementation of FFTW3. Proceedings of the IEEE, 93(2), 216–231. https://doi.org/10.1109/JPROC.2004.840301
11. Frigo, M., Leiserson, C. E., Prokop, H., & Ramachandran, S. (1999). Cache-oblivious algorithms. In 40th Annual Symposium on Foundations of Computer Science (pp. 285–297). IEEE. https://doi.org/10.1109/SFFCS.1999.814600
12. Hairer, E., & Wanner, G. (1996). Solving ordinary differential equations II: Stiff and differential-algebraic problems (2nd ed.). Springer. https://doi.org/10.1007/978-3-642-05221-7
13. Kapoor, S., & Narayanan, A. (2023). Leakage and the reproducibility crisis in machine-learning-based science. Patterns, 4(9), Article 100804. https://doi.org/10.1016/j.patter.2023.100804
14. Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., & Yang, L. (2021). Physics-informed machine learning. Nature Reviews Physics, 3(6), 422–440. https://doi.org/10.1038/s42254-021-00314-5
15. Kaufman, S., Rosset, S., Perlich, C., & Stitelman, O. (2012). Leakage in data mining: Formulation, detection, and avoidance. ACM Transactions on Knowledge Discovery from Data, 6(4), Article 15. https://doi.org/10.1145/2382577.2382579
16. Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. Annual Review of Earth and Planetary Sciences, 26, 643–696. https://doi.org/10.1146/annurev.earth.26.1.643
17. Mignan, A., & Broccardo, M. (2019). One neuron versus deep learning in aftershock prediction. Nature, 574(7776), E1–E3. https://doi.org/10.1038/s41586-019-1582-8
18. Mignan, A., & Broccardo, M. (2020). Neural network applications in earthquake prediction (1994–2019): Meta-analytic and statistical insights on their limitations. Seismological Research Letters, 91(4), 2330–2342. https://doi.org/10.1785/0220200021
19. Mousavi, S. M., & Beroza, G. C. (2022). Deep-learning seismology. Science, 377(6607), Article eabm4470. https://doi.org/10.1126/science.abm4470
20. Mousavi, S. M., Sheng, Y., Zhu, W., & Beroza, G. C. (2019). STanford EArthquake Dataset (STEAD): A global data set of seismic signals for AI. IEEE Access, 7, 179464–179476. https://doi.org/10.1109/ACCESS.2019.2947848
21. Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Köpf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., … Chintala, S. (2019). PyTorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems 32 (pp. 8024–8035). https://arxiv.org/abs/1912.01703
22. Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 686–707. https://doi.org/10.1016/j.jcp.2018.10.045
23. Rice, J. R., & Ruina, A. L. (1983). Stability of steady frictional slipping. Journal of Applied Mechanics, 50(2), 343–349. https://doi.org/10.1115/1.3167042
24. Rouet-Leduc, B., Hulbert, C., Lubbers, N., Barros, K., Humphreys, C. J., & Johnson, P. A. (2017). Machine learning predicts laboratory earthquakes. Geophysical Research Letters, 44(18), 9276–9282. https://doi.org/10.1002/2017GL074677
25. Rucker, C., & Erickson, B. A. (2024). Physics-informed deep learning of rate-and-state fault friction. Computer Methods in Applied Mechanics and Engineering, 430, Article 117211. https://doi.org/10.1016/j.cma.2024.117211
26. Ruina, A. (1983). Slip instability and state variable friction laws. Journal of Geophysical Research: Solid Earth, 88(B12), 10359–10370. https://doi.org/10.1029/JB088iB12p10359
27. Scholz, C. H. (1998). Earthquakes and friction laws. Nature, 391(6662), 37–42. https://doi.org/10.1038/34097
28. SEG Technical Standards Committee. (2017). SEG-Y_r2.0: SEG-Y revision 2.0 data exchange format (R. Hagelund & S. A. Levin, Eds.). Society of Exploration Geophysicists. https://seg.org/wp-content/uploads/2025/11/seg_y_rev2_0_mar2017.pdf
29. Singh, A. (2026). torch-friction (Version 0.1.0) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.18366445
30. Trnkoczy, A. (2012). Understanding and parameter setting of STA/LTA trigger algorithm. In P. Bormann (Ed.), New manual of seismological observatory practice 2 (NMSOP-2). Deutsches GeoForschungsZentrum GFZ. https://doi.org/10.2312/GFZ.NMSOP-2_IS_8.1
31. Withers, M., Aster, R., Young, C., Beiriger, J., Harris, M., Moore, S., & Trujillo, J. (1998). A comparison of select trigger algorithms for automated global seismic phase and event detection. Bulletin of the Seismological Society of America, 88(1), 95–106. https://doi.org/10.1785/BSSA0880010095
32. Yilmaz, Ö. (2001). Seismic data analysis: Processing, inversion, and interpretation of seismic data (2nd ed.). Society of Exploration Geophysicists. https://doi.org/10.1190/1.9781560801580

Code & archives

• torch-friction (Zenodo DOI)

• Basalt (Zenodo DOI)

• fast-trigger

• seismic-linter (PyPI)

This article was first published at amansinghnp.com.