gravimetry

Turn monthly gravity fields from GRACE and GRACE-FO into maps of water storage, and model the gravity change an earthquake causes. Built on geopotential.

Water stored on land and in the ocean in July 2026, against the 2004–2009 average.

Install

cargo add gravimetry geopotential

Example: from a monthly field to water height

The standard chain of the GRACE Tellus maps (Landerer & Swenson 2012):

1. subtract a baseline;
2. remove the stripes;
3. smooth over 300 km;
4. convert to centimetres of water.

use geopotential::{Coefficients, Functional};

let (gm, radius) = (3.986004415e14, 6_378_136.3);

// A month and the baseline it is compared with (here: one changed coefficient).
let mut month = Coefficients::zeros(gm, radius, 60);
month.set(2, 0, 1e-10, 0.0);
let baseline = Coefficients::zeros(gm, radius, 60);

let change = gravimetry::anomaly(&month, &baseline);
let change = gravimetry::destripe(&change);
let change = gravimetry::smooth(&change, &gravimetry::gaussian_weights(300e3, radius, 60));

let love = gravimetry::love_numbers(60).unwrap();
let water = change.grid(&Functional::WaterHeight { love: &love }, radius, &[90f64.to_radians()], &[0.0]);
println!("{:.0} mm of water at the North Pole", water[0] * 1000.0); // 19 mm

With real data, read the month and the technical notes with gravity-field. The full chain also replaces C̄20, C̄30 and degree 1.

Example: the gravity change of an earthquake

Left: the gravity step the July 2025 Kamchatka earthquake left in GRACE-FO. Right: the USGS fault model through the same filter, computed with this crate.

use gravimetry::coseismic::{Fault, Medium};

// A thrust fault: 100 × 50 km, 8 m of slip, its centre 20 km deep.
let fault = Fault {
    lat: 52.5, lon: 160.0, depth: 20e3,
    strike: 217.0, dip: 15.0, rake: 90.0,
    length: 100e3, width: 50e3, slip: 8.0, opening: 0.0,
};

// Uplift and gravity change at the surface straight above it.
let above = fault.surface_change(0.0, 0.0, &Medium::default());
println!("uplift {:.2} m, gravity {:+.0} µGal", above.uplift, above.gravity * 1e8); // 1.20 m, +115 µGal

// What a satellite sees: surface masses on a grid, as coefficients, smoothed like the data.
let lats: Vec<f64> = (0..50).map(|i| (47.0 + 0.2 * i as f64).to_radians()).collect();
let lons: Vec<f64> = (0..70).map(|j| (153.0 + 0.2 * j as f64).to_radians()).collect();
let at_sea = |_lat: f64, _lon: f64| true;
let masses = gravimetry::coseismic::surface_masses(&[fault], &lats, &lons, &Medium::default(), at_sea);
let field = gravimetry::mass_coefficients(&masses, 3.986004415e14, 6_378_136.3, 60);
let field = gravimetry::smooth(&field, &gravimetry::gaussian_weights(300e3, field.radius, 60));

For a real earthquake, coseismic::read_coulomb reads the coulomb.inp file that the USGS publishes with every finite-fault model.

What is inside

How it is checked

• Love numbers: against the published table.
• Smoothing weights: against their closed-form limit.
• Stripe removal: on synthetic patterns.
• Fault model:
  • Okada's (1985) published check values;
  • mass conservation: a buried fault's equivalent mass sums to zero;
  • on the 2025 Kamchatka fault model, its uplift matches the USGS's own surface displacement with r = 0.89.
• Point masses: their potential matches Newton's law.

Not included yet: glacial isostatic adjustment, DDK filters, basin averages and leakage corrections.

Part of quantum-gravimeter.