Crate gravimetry 

gravimetry

Satellite gravimetry from instrument to map. Model the satellites, from today’s GRACE-FO to missions with quantum (cold-atom) sensors, and simulate what gravity field they recover. Then turn the monthly fields of 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 each month with geopotential::read_gfc and the low-degree corrections with technotes::read_tn14 and read_tn13. The full chain replaces C̄20, C̄30 and degree 1 before subtracting the baseline.

Example: a quantum accelerometer against GRACE-FO

Instrument noise, reproducing Figure 1 of Zingerle et al. (2024) from its equations.

use gravimetry::atom::AtomInterferometer;
use gravimetry::missions::{Accelerometer, Ranging, pair_noise};

// A cold-atom accelerometer like CARIOQA's: 1 s between pulses, 80% contrast,
// 100,000 atoms per shot, one shot every 3.5 s, double diffraction in orbit.
let atoms = AtomInterferometer { t: 1.0, contrast: 0.8, atoms: 1e5, cycle: 3.5, diffraction: 2 };
let white = atoms.acceleration_asd(); // about 2.3e-10 m/s² per √Hz

// Noise of a satellite pair at 0.1 mHz, the longest wavelengths, in m/s² per √Hz.
let quantum = pair_noise(Accelerometer::AtomInterferometer { white }, Ranging::Lri2033, 1e-4);
let today = pair_noise(Accelerometer::SuperStar, Ranging::GraceFo, 1e-4);
println!("quantum {quantum:.1e}, GRACE-FO {today:.1e}"); // 3.2e-10 against 1.0e-9
assert!(quantum < today / 3.0);

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.

What is inside

Item	What it does
missions	Noise spectra of electrostatic and cold-atom accelerometers and laser ranging, and of a whole satellite pair (Zingerle et al. 2024, eqs. 1–8)
atom	Phase, shot noise, and accelerometer and gradiometer noise of a cold-atom sensor, in orbit (double diffraction) or on the ground (single), its phase sensitivity and the gravitational Aharonov–Bohm rate (Lévèque et al. 2021; Stray et al. 2022; CARIOQA; Trimeche et al. 2019)
gradiometry	The error of each coefficient a gradiometer mission like GOCE achieves
local	Gravity and gravity gradient of a buried sphere and horizontal cylinder, the shapes a ground survey fits, and the averaging time an instrument needs to see them
clocks	Clocks as height meters: frequency ↔ geopotential ↔ height, the clock rates of TT and TCG, proper time along an orbit, and the GNSS relativistic corrections (eccentricity, Sagnac, Shapiro) (IERS Conventions 2010, ch. 10; IS-GPS-200)
tides	The lunar and solar tide at a gravimeter, the gravimetric and diminishing factors of the Love numbers, and the permanent tide that separates mean-tide, zero-tide and tide-free gravity and geoids (Longman 1959, as LongmanTide codes it; Mäkinen 2021)
sim	A closed-loop simulation of a GRACE-type pair: simulate, estimate the field back, and compare its errors with the predicted ones
noise	Noise with any spectrum, its whitening, and the FFT
technotes	NASA’s TN-14 (C̄20 and C̄30 from laser ranging, Loomis et al. 2020) and TN-13 (degree 1, Sun et al. 2016), from PO.DAAC
love_numbers	How much the solid Earth sags under a load, per degree (Wahr et al. 1998, Table 1)
gaussian_weights, smooth	Averaging over a chosen radius (Wahr et al. 1998, after Jekeli 1981)
destripe	Removes the north–south stripes of GRACE noise (Swenson & Wahr 2006), ported from gravity-toolkit
mean, anomaly	Baseline means and differences, across GM and radius conventions
mass_coefficients	Spherical-harmonic coefficients of point masses
transport_weights	The smallest mass flow that carries a change between months
coseismic	Uplift and gravity change of faults in an elastic half-space (Okada 1985; Okubo 1992); reads USGS coulomb.inp files and corrects for seawater pushed aside by seafloor uplift
kalman	A Kalman smoother for series with gaps and per-month errors, and a maximum-likelihood fit of its settings (Rauch et al. 1965)
atmosphere	Sound, buoyancy and Lamb waves in the air, the mass of a pressure change, and the ocean’s forced response to a fast pressure pulse (Wright et al. 2022; Proudman 1929)

The simulated error of four mission classes, from today’s to a quantum pair, against the errors of a real month.

Applications

• Drought and groundwater. GRACE and GRACE-FO mass change feeds weekly drought indicators and statutory groundwater reporting; the processing chain here is the one those products use.
• Ice sheets and sea level. examples/ice_and_sea_level.rs turns Greenland’s published mass trend into its sea-level contribution.
• Earthquakes. coseismic models the gravity change of a fault and reads USGS finite-fault files; prompt_lead gives how far ahead of the P-wave a light-speed gravity signal arrives.
• Site surveys and monitoring. local gives the gravity and gradient of a buried void or dense body, and the averaging time a quantum gravimeter needs to resolve it.
• Absolute gravimetry. tides corrects a reading for the lunar and solar tide and moves it between the zero-tide convention of gravity and the mean-tide convention of the IHRS.
• Chronometric levelling. clocks turns an optical-clock comparison into a height difference, and gives a GNSS receiver its relativistic corrections.
• Mission design. missions, atom and sim size the sensors and simulate what field a satellite pair recovers, including cold-atom instruments.

How it is checked

• Noise models: they reproduce Figure 1 of Zingerle et al. (2024).
• Cold-atom sensors: they reproduce the published CARIOQA, GRICE and Trimeche et al. (2019) sensitivities.
• Gradiometry errors: the closed form matches brute-force least squares.
• Simulator: it recovers a known field exactly without noise; with noise, a χ² test confirms that its actual errors match its predicted ones.
• Technical notes: real NASA files, excerpted in tests/data.
• 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.
• Clocks: the Tokyo Skytree (452.6 m is 21.18 Hz on strontium, and GRS80 normal gravity there is the measured g to 2×10⁻⁶); the JILA and Wisconsin redshift gradients; the Sun’s 633.10 m/s; Cassini’s 262.4 µs Shapiro delay; the GPS clock offset 4.4647×10⁻¹⁰ and IS-GPS-200’s F.
• Tides: LongmanTide’s test vector to 3×10⁻¹¹ mGal; averaged over the 18.61-year nodal cycle, Longman’s tide is δ₂ times the permanent tide to 0.2 µGal; the permanent tide matches Mäkinen’s (2021) Eqs. 19–21 to their printed digits.
• Kalman smoother: it matches the exact Gaussian-process solution to 1e-10, its fit recovers the settings of simulated data, and 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), with standard gravity; 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; time-variable gravity in the simulator, so no temporal aliasing.

One of eight crates in physica: physica (all of them, plus cross-domain analogies) · metrology (constants) · geopotential (the gravity field) · gravimetry (missions and GRACE) · gravitational (waves) · gravitics (frontier gravity) · magnetics (magnetism and MHD) · plasmoid (fusion plasma).

Modules

atmosphere

Waves in the atmosphere that a gravity satellite can feel: sound, buoyancy (“gravity”) waves and the Lamb wave, and the mass a pressure change puts above a square metre of ground.

atom

Cold-atom (light-pulse) interferometers as space accelerometers and gradiometers, at the quantum-projection-noise limit.

clocks

Clocks as gravimeters. A clock higher in the Earth’s field runs faster by gΔh/c², 1.1×10⁻¹⁶ per metre, so comparing two clocks measures their geopotential difference, and an optical clock at 10⁻¹⁸ resolves a centimetre of height (chronometric levelling). The same physics sets the rates of the relativistic time scales and the corrections a GNSS receiver applies.

coseismic

Gravity change from an earthquake: rectangular dislocations in an elastic half-space (Okubo 1992), the USGS coulomb.inp finite-fault format, and the equivalent surface masses a satellite sees.

gradiometry

From instrument noise to gravity-field errors: formal errors of spherical-harmonic coefficients for a gravity-gradiometer mission, and the spatial resolution of a degree.

kalman

A Kalman smoother for monthly series with gaps, such as one gravity coefficient over the GRACE and GRACE-FO record.

local

Forward models for a buried body: the gravity and gravity gradient of a sphere and of a horizontal cylinder, the shapes a ground survey actually fits.

missions

Gravity-mission sensors and their noise: accelerometers, inter-satellite ranging, and the along-track noise of a satellite pair.

noise

Time series with a given noise spectrum, and whitening of series by that spectrum: the two frequency-domain steps of a closed-loop simulation.

sim

A minimal closed-loop simulation of a GRACE-type satellite pair.

technotes

GRACE / GRACE-FO technical-note replacement series: TN-14 (C̄20, C̄30 from satellite laser ranging) and TN-13 (degree-1 / geocentre coefficients). Both are anonymous downloads from the PO.DAAC docs bucket. Epochs are Julian years, as for read_gfc.

tides

The Earth tide at a gravimeter: Longman’s lunar and solar tide, the Love-number factors that scale a rigid Earth’s tide to the observed one, and the permanent tide that separates the mean-tide, zero-tide and tide-free conventions.

Constants

OCEAN_AREA

Area of the world ocean, m²: 3.618×10⁸ km², so 361.8 Gt of water raise sea level 1 mm.

Functions

anomaly

field − reference, with the reference rescaled to the field’s GM and radius, to the lower of the two degrees.

destripe

Remove north–south stripes (correlated errors) the way Swenson & Wahr (2006, GRL 33, L08402) do: for each order m ≥ 5, fit a moving quadratic across degree, separately for even and odd degrees, over 2nₛ+1 coefficients with nₛ = max(round(15·e^(−m/10)), 2), and subtract the fit. Sequences too short for one window are removed entirely. Apply it to an anomaly (typically truncated to degree 60), before Gaussian smoothing.

gaussian_weights

Gaussian smoothing weights for averaging radius radius (m; the distance at which the spatial kernel falls to half), on a sphere of radius earth_radius, for n = 0..=nmax. Jekeli (1981) / Wahr et al. (1998, eqs. 32–34), normalized so W_0 = 1.

love_numbers

Elastic load Love numbers k'_n (PREM, Han & Wahr 1995) as tabulated by Wahr et al. (1998, JGR 103, Table 1), linearly interpolated between tabulated degrees (error < 0.05%), for n = 0..=nmax. k'_1 = 0.027 is the centre-of-figure value. None beyond degree 200, the end of the table.

mass_coefficients

Spherical-harmonic coefficients, to degree nmax, of point masses (lat, lon, kg) (radians) on the sphere of radius radius, in the convention of a field with gm: C̄nm + iS̄nm = (m_k / M) · P̄nm(sin φ_k) · e^{imλ_k} / (2n + 1) summed over the masses, with M = gm / G. Degree 0 holds the net mass; it is kept so the caller can check it.

mean

The mean of several fields, in the first field’s GM and radius (others are rescaled), to the lowest degree among them. None if fields is empty.

sea_level_equivalent

Global mean sea-level rise (m) from mass (kg) of water added to the ocean and spread evenly: M/(ρ_w A). Ignores self-attraction and loading (sea-level fingerprints).

smooth

Multiply every degree n by weights[n] (e.g. gaussian_weights).

transport_weights

Degree weights R²/(n(n+1)) (zero at degree 0) that turn the rate of change of a field into its minimum-norm mass-transport potential φ, the solution of ∇²φ = −dσ/dt on the sphere. Apply them with smooth to a field of yearly change, synthesize geopotential::Functional::WaterHeight and multiply by the water density: that is φ in kg/yr, and its horizontal gradient is the smallest mass flux (kg per metre per year) that could carry the change. It points toward where mass is gathering.