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.

Example: fill the gaps in a monthly record

The GRACE record has gaps: missing months, and almost a year between the two missions. The Kalman smoother treats each coefficient as a slowly wandering signal plus noise. It gives the best estimate and its uncertainty for every month, gaps included.

use gravimetry::kalman::{Process, smooth};

// Five months, the middle one missing, each observed with an error of 1.
let times = [0.0, 1.0, 2.0, 3.0, 4.0].map(|month| month / 12.0);
let observed = [Some((1.0, 1.0)), Some((2.0, 1.0)), None, Some((2.5, 1.0)), Some((2.0, 1.0))];

// A signal that varies by about 2 over about 6 months.
let result = smooth(&times, &observed, Process { sd: 2.0, tau: 0.5 });
println!("month 3: {:.2} ± {:.2}", result.mean[2], result.sd[2]);

kalman::fit picks the signal's size and correlation time, and checks the stated errors, by maximum likelihood over many series at once.

The smoother fitted to 24 years of monthly fields: how far COST-G's stated errors need scaling per degree, and how well withheld months are predicted.

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.
• Kalman smoother:
  • it matches the exact Gaussian-process solution to 1e-10;
  • its fit recovers the settings of simulated data;
  • on the real record, withheld months are predicted with calibrated errors.
• Atmosphere: the isothermal values at 250 K (sound 317 m/s, cut-off 0.0217 rad/s, buoyancy 0.0196 rad/s); the Lamb wave is trapped at every wavelength; the Tonga Lamb speed of 318.2 m/s is sound at 252 K; a slow pressure pulse gives the inverse barometer.

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

Part of quantum-gravimeter.