Seafloor compliance, pressure-gauge calibration and Bayesian shear-velocity inversion for broadband ocean-bottom seismometers.
ComPy turns raw ocean-bottom seismometer (OBS) and differential pressure gauge (DPG) records into crustal shear-velocity models. It covers the full chain: data cleaning, tilt and coherent-noise removal, Rayleigh-wave pressure calibration, stacked compliance estimation and Bayesian inversion for Vs and serpentinization. It was developed for the RHUM-RUM experiment in the Indian Ocean (Aminian et al. 2025, GJI) and runs unattended on any broadband OBS network stored as an SDS archive.
- Highlights
- Method
- Installation
- Quick start
- Workflow, step by step
- Results
- Repository structure
- Documentation
- Citation
- License and acknowledgments
- End-to-end pipeline. One command takes an SDS archive and StationXML to
compliance, Vs models, figures and a report (
pipeline/run_pipeline.py). - Stacked compliance with full error budget. Cross-spectra stacked over non-overlapping windows, gravity and wave-attraction corrections, Bendat & Piersol errors and monthly stacks for the covariance.
- Rayleigh-wave DPG calibration. Per-event gains from teleseismic Rayleigh
waves, cross-checked against microseism and local-P estimates
(
compy_pressure). - Vs inversion that says what it resolves. One command (
compy-invert) from a compliance curve to Vs with credible intervals: errors measured from the station's own windows, a prior about a reference model, Vp, density, Moho depth and pressure gain uncertain, sampling repeated until it converges, and each depth marked as resolved by the data or not (docs/INVERSION.md). - Two networks processed. RHUM-RUM (YV, 8 stations, Indian Ocean) and Pacific Array Young ORCA (XE, 30 stations, Pacific).
- Tested. CI on Python 3.10–3.13 with lint, format, wheel build and an offline multi-chain benchmark.
Long-period seafloor pressure from infragravity waves deforms the crust. The normalized compliance is the vertical displacement per unit pressure, scaled by the water-wave wavenumber:
where
Python 3.10 or later. With conda (recommended, for Numba):
git clone https://github.com/MohammadAmin-Aminian/ComPy.git
cd ComPy
conda env create -f environment.yml
conda activate compyor with pip:
python -m venv .venv && source .venv/bin/activate
python -m pip install -e '.[dev,pressure,fast]'
python -m pytest -qThe distribution is seafloor-compy; it installs the compy package and
the compy_pressure package. The ComPy 2.x module names (inv_compy,
ffplot, compy_processing, ...) still import, with a deprecation warning
that names the new location. Dependencies: NumPy, SciPy, Matplotlib, ObsPy,
TiSKitPy 2.3 and disba; Numba is optional and speeds up the inversion about
140×.
Offline synthetic inversion (seconds):
MPLBACKEND=Agg python examples/synthetic_inversion.pyA full network run from an SDS archive:
export COMPY_RUN_DIR=~/runs/xe # metadata/, work/, results/, logs/
python pipeline/run_pipeline.py CC11 CC05 --workers 8 --processes 8 --fetchThe pipeline runs stages A → G (preprocessing, cleaning, calibration,
compliance, inversion, figures, report) with checkpoints and a JSON run
manifest. All network-specific values live in pipeline/config.py; see the
pipeline guide.
The figures follow one station through the whole pipeline: XE.CC11 (Pacific Array Young ORCA, 7.54°S 133.42°W, 4636 m water depth, April 2018 – May 2019). Steps 1 and 8 place it within the XE network.
| Step | Stage | What it does |
|---|---|---|
| 1. Network and data | A | Read waveforms and StationXML responses |
| 2. Response removal and cleaning | A, B | Decimate, remove responses, drop clipped data, mask earthquakes |
| 3. Tilt and coherent noise | B | Remove tilt and horizontal-coherent noise from Z |
| 4. Pressure-gauge calibration | D | Rayleigh-wave gain of the DPG |
| 5. Coherence and window selection | C | Select windows with high Z–P coherence |
| 6. Compliance | C | Stacked compliance with errors and monthly stacks |
| 7. Bayesian Vs inversion | E | Multi-chain inversion for Vs, Vp/Vs and gain |
| 8. Vs profiles and serpentinization | F | Station values and network profiles |
Waveforms are read from an SDS archive with their full StationXML responses; response epochs, units, timing and gaps are checked per station. XE has 30 Scripps broadband OBS (Trillium 240 + DPG, 50 Hz) at 4386–4824 m water depth on ~40 Ma Pacific seafloor. CC11 is in the central (CC) line.
python pipeline/fetch_inventory.py && python pipeline/fetch_catalog.py
python pipeline/a_preprocess.py CC11 --workers 8
XE network. Red: inverted (20 stations, including CC11); orange: compliance only; grey: not processed.
Each channel is decimated with an anti-alias filter and the instrument response is removed with an explicit pre-filter: acceleration or displacement for the seismometer, pascals for the DPG. Ten-minute blocks containing a clipped or flat sample are removed from all channels. Earthquakes are masked from the global catalogue with TiSKitPy, and STA/LTA picks local events. No hourly glitches were found on XE (fold strength below 2.0), so glitch removal is off for this network.
import tiskitpy
spans = tiskitpy.TimeSpans.from_eqs(
(z.stats.starttime, z.stats.endtime),
minmag=5.5,
days_per_magnitude=0.5,
save_eq_file=False,
)
XE.CC11 spectra. (a) Vertical acceleration after response removal (blue) and after cleaning (orange), median and 10–90 %, with the Peterson noise models. (b) Calibrated pressure. Green: compliance band.
Hourly tilt angles from TiSKitPy's CleanRotator are median-smoothed and
interpolated per sample. Horizontal-coherent noise is then removed from the
vertical only, in 24-h windows that avoid the earthquake spans.
import compy
from compy import cleaning
rotated = cleaning.rotate_smoothly(st, hour_centres, angles, azimuths)
local = cleaning.sta_lta_spans(rotated.select(component="Z")[0])
cleaned = compy.remove_coherent_noise(rotated, 24, avoid_spans=spans, cover_all=True)
XE.CC11 tilt. Hourly azimuth and incidence angle that minimize vertical variance, coloured by variance reduction.
The DPG gain is fitted per teleseism by comparing the measured pressure-to-acceleration ratio of the Rayleigh wave with the water-layer prediction, in bins of high coherence. The station gain is the median of the accepted events. For CC11 it is 0.698 ± 0.012 from 12 events, against 0.702 from the independent BRUIT-FM microseism method.
python examples/3_DPG_Calibration.py SDS/ stationxml/ XE CC11 events.geojson
XE.CC11 calibration. (a) P–Z coherence of 12 teleseisms (M ≥ 7). (b) Measured, theoretical and corrected spectral ratios; gain g = 0.698 ± 0.012.
The cleaned record is cut into non-overlapping 35-min windows. A window is kept when its median Z–P coherence in the compliance band exceeds 0.9 and its pressure and acceleration levels lie inside the selection gates. CC11 keeps 7314 windows over 14 months.
XE.CC11, February 2019. Pressure and vertical spectrograms and the Z–P coherogram; dashed lines bound the compliance band.
Compliance comes from cross-spectra stacked over all selected windows, with the gravity and wave-attraction corrections and the Bendat & Piersol error. Monthly stacks give the data covariance and track changes in time. For CC11 the band is 4.9–14.6 mHz, the mean coherence 0.97, and the monthly stacks stay within −0.9 % and +1.2 % of the full stack.
from compy.compliance import GJI2025, stacked_compliance
pair = cleaned.select(channel="BHZ") + cleaned.select(channel="BDH")
windows = compy.split_stream(pair, 2100) # 35-min windows
result = stacked_compliance(
windows, 4636, 0.698, band=(0.0049, 0.0146), criteria=GJI2025, nseg=1024
)
result.compliance, result.uncertainty, result.n_windows
XE.CC11 selected windows (150 of 7314 shown): corrected vertical acceleration, calibrated pressure, coherence and per-window compliance with the stacked estimate (red).
Since ComPy 3 the default is compy.inversion.auto (docs/INVERSION.md;
tests in docs/inversion_overhaul_2026-11/REPORT.md).
Compliance measures stiffness, so Vs follows only with Vp and density: both are
tied to Vs by rock relations with sampled factors. A soft-sediment unit, crustal
layers that double in thickness, a free Moho depth and a mantle half-space are
sampled with a prior about the CRUST1.0 model of the station (or a model you
give). Errors come from the bootstrap of the station's own windows plus terms for
the slope and level of the curve. Eight chains with parallel tempering run until
they converge.
compy-invert CC11_compliance_compy.npz --lat -7.5368 --lon -133.4166 --crust1 ~/crust1 --fetch-crust1from compy.inversion import auto, dram
model = auto.model_from_crust1(dram.crust1_column(-7.5368, -133.4166, "crust1"))
result = auto.invert(data, model, gain_sigma=0.05) # data: auto.compliance_data(...)
result.quality(), result.profile() # label with reasons; Vs, prior, resolved depths
XE.CC11 with ComPy 3. Data and fit, residuals in units of the error, Vs with its prior and the reference, and how much the data narrowed each parameter. The black bar marks the depths the data resolve.
The numbers and figures of this section and of the Results below were made with ComPy 2.4 to 2.8. ComPy 3 changes Vs (by about +6 % at 2 km and +11 % at 5 km in the median over 155 stations, more at the YV stations) and the sediment thickness; see docs/inversion_overhaul_2026-11/REPORT.md. They are kept as the record of those runs.
Posterior draws are converted to serpentinization with the Miller & Christensen (1997) relation used in Aminian et al. (2025, Fig. 10c): 0 % at 4.2 km/s and 100 % at 2.2 km/s. The conversion attributes all of the Vs reduction to serpentinization, so it gives an upper bound.
beta_q = post.serpentinization_quantiles(z) # (5 quantiles, len(z)) in %XE.CC11, median (5–95 %):
| Depth | Vs [m/s] | Serpentinization upper bound [%] |
|---|---|---|
| 0.5 km | 2530 (1942–3428) | 84 (39–100) |
| 1 km | 2848 (2387–3369) | 68 (42–91) |
| 2 km | 3168 (2890–3482) | 52 (36–66) |
| 4 km | 3496 (3261–3976) | 35 (11–47) |
Sediment thickness 121 m (20–182), crustal Vp/Vs 2.17 (1.88–2.37).
XE network. (a) Median Vs by sub-array (CC11 in Central). (b) Pooled posterior with free Vp/Vs and with Vp/Vs = 1.85. (c) Serpentinization upper bound.
30 stations on ~40 Ma seafloor, 20 inverted, 18 with a Rayleigh-wave gain.
- Compliance at 4.9–14.6 mHz with Z–P coherence 0.93–0.97 and thousands of windows per station; month-to-month variation of 1–3 % (CC07: ±10–25 %).
- 80–170 m of sediment, then Vs of about 2.4 km/s at 0.3 km, 2.9 km/s at 1 km and 3.2 km/s at 2 km.
- Rayleigh-wave gains agree with the independent BRUIT-FM microseism gains (median ratio ≈ 0.99).
XE compliance by sub-array, ±1σ of the full covariance.
Full report (PDF, LaTeX source and tables): docs/results/xe_young_orca_2026-10-29.
Eight broadband stations, 2012–2013, the network of Aminian et al. (2025). Final inversion: free Vp/Vs, gain sampled, stable and post-jump epochs inverted together. Vs in m/s, median (5–95 %).
| Station | Setting | χ²/N | R-hat | Sediment [m] | Vs 0.5 km | Vs 1 km | Vs 2 km |
|---|---|---|---|---|---|---|---|
| RR28 | Somali plate | 2.9/6 | 1.012 | 311 (138–398) | 2157 (1450–3760) | 2789 (2150–3725) | 3421 (2910–3931) |
| RR29 | Somali plate | 2.9/5 | 1.005 | 91 (15–139) | 3186 (2333–4252) | 3426 (2709–4217) | 3678 (3106–4218) |
| RR34 | Somali plate | 2.4/6 | 1.008 | 219 (50–309) | 2258 (1723–3727) | 2811 (2203–3560) | 3185 (2757–3686) |
| RR36 | SWIR | 6.8/16 | 1.017 | 72 (10–143) | 2236 (1756–3320) | 2578 (2082–3094) | 3090 (2740–3594) |
| RR38 | SWIR | 33.0/12 | 1.015 | 224 (138–251) | 3062 (2307–4369) | 3722 (3025–4402) | 4058 (3697–4471) |
| RR40 | SWIR | 4.2/10 | 1.009 | 164 (40–219) | 2964 (2092–4165) | 3269 (2581–3978) | 3556 (3078–3994) |
| RR50 | Rodrigues TJ | 5.7/12 | 1.005 | 80 (14–112) | 3389 (2502–4346) | 3672 (2980–4363) | 3964 (3495–4429) |
| RR52 | CIR | 3.3/9 | 1.007 | 46 (6–81) | 2823 (2332–3632) | 3193 (2779–3672) | 3590 (3283–4013) |
RR38's pressure channel changes behaviour after its sensitivity jump, so RR38 is excluded from the interpretation. Figures, variants and full tables: docs/results/yv_inversion_2026-10-29 and docs/results/yv_compliance_2026-10-28.
ComPy/
├── compy/ # the package
│ ├── streams.py # window splitting and stream checks
│ ├── cleaning.py # smooth tilt rotation, periodic transients, STA/LTA spans
│ ├── preprocessing.py # tilt correction, coherent-noise removal
│ ├── compliance.py # spectra, window selection, stacked compliance, errors
│ ├── inversion/ # Vs inversion
│ │ ├── auto.py # default inversion (ComPy 3): invert, SeafloorModel, compy-invert
│ │ ├── errors.py # covariance of a stacked compliance curve from its windows
│ │ ├── dram.py # DRAM sampler, forward model, CRUST1.0, the 2.x model, serpentinization
│ │ ├── metropolis.py # Metropolis sampler of the GJI paper, tempering
│ │ ├── occam.py # Occam smooth-model inversion
│ │ ├── minimum_structure.py # minimum-structure inversion after Crawford (2004)
│ │ ├── literature_mcmc.py # Monte Carlo with Zha & Webb (2016), Zheng et al. (2024)
│ │ ├── impedance.py # independent forward model (surface impedance)
│ │ ├── propagator.py # forward model (Crawford's raydep) and kernels
│ │ └── models.py # starting and reference models
│ ├── diagnostics.py # R-hat, ESS, reproducible run records
│ ├── plotting/ # style, spectra, inversion and QC figures
│ └── numerics.py # constants, gravity-wave dispersion, shared maths
├── inv_compy.py, ffplot.py, compy_*.py # ComPy 2.x names (deprecated aliases)
├── compy_pressure/ # DPG QC, Rayleigh/microseism/local-P calibration, noise models
├── pipeline/ # end-to-end network pipeline (stages A–G)
├── examples/ # offline synthetic examples and per-step scripts
├── tests/ # pytest suite
├── docs/
│ ├── USER_GUIDE.md # API reference
│ ├── INVERSION.md, PRESSURE.md, VALIDATION.md, ...
│ ├── images/ # logos and figures used by the documentation
│ ├── benchmark/ # CI inversion benchmark record
│ └── results/ # processed networks: scripts, figures, tables
│ ├── xe_young_orca_2026-10-29/
│ ├── yv_compliance_2026-10-28/
│ └── yv_inversion_2026-10-29/
└── pressure/ # pressure-gauge campaign scripts, method notes and tests
| Document | Content |
|---|---|
| Documentation index | All guides, method notes and results |
| User guide | Input conventions, estimators, inversion API, pressure QC, plotting |
| Pipeline guide | Running a whole network, stages and configuration |
| Examples | Offline synthetic examples and per-step scripts |
| Inversion notes | Parameterization, sampler and validation |
| Alternative inversions | The Moho rule, published Monte Carlo setups, Crawford's minimum-structure inversion |
| Pressure tools | compy_pressure commands and outputs |
| Validation | Test coverage and scientific references |
| Changelog | Release history |
If you use ComPy, please cite the article and the software version you used (CITATION.cff, Zenodo).
@article{aminian2025shallow,
author = {Aminian, Mohammad Amin and Crawford, Wayne and Stutzmann, Éléonore
and Montagner, Jean-Paul and Cannat, Mathilde and Hadziioannou, Céline},
title = {Shallow crustal structures of the Indian ocean derived from compliance function analysis},
journal = {Geophysical Journal International},
year = {2025},
volume = {242},
number = {3},
pages = {ggaf253},
doi = {10.1093/gji/ggaf253}
}ComPy is released under the GPL-3.0. Contributions are welcome; see CONTRIBUTING.md.
Author: Mohammad-Amin Aminian · LinkedIn · Research software portfolio
Developed at the Institut de Physique du Globe de Paris with support from the SPIN Innovative Training Network, funded by the European Commission under Horizon 2020 Marie Skłodowska-Curie Actions.
- OBS Transient Cleaner: periodic OBS instrument-transient removal.
- ComPy Inversion Tuner: tuning of compliance-inversion controls.
- RHUM-RUM Geospatial Mapper: bathymetry, OBS-network and tectonic maps.
- VRE Seismic Enhancement: Virtual Resolution Enhancement for seismic sections.
- Gabor Seismic Filter: orientation-selective 2-D seismic filtering.