Gravity Remote Sensing (GRACE)
modelling Level 4

Gravity Remote Sensing (GRACE)

How do we measure water storage changes and ice loss from space? GRACE (Gravity Recovery and Climate Experiment) satellites detect tiny gravity variations caused by mass redistribution. This model derives gravitational potential equations, implements mass anomaly calculations, demonstrates spherical harmonic analysis, and shows applications from groundwater depletion to ice sheet mass balance.

Prerequisites: gravitational potential, mass anomaly, spherical harmonics, water equivalent thickness

Updated 10 min read

1. The Question

How much water is California losing during drought?

GRACE mission (2002-2017, GRACE-FO 2018-present):

Twin satellites measure gravity variations.

Principle:

Mass changes → gravity changes → satellite separation changes

Measurement:

Inter-satellite distance via microwave ranging.

Precision: ~10 μm range change → ~1-2 cm equivalent water thickness

Temporal resolution: Monthly

Spatial resolution: ~300-400 km

Applications:

  • Ice sheet mass balance (Greenland, Antarctica)
  • Groundwater storage changes
  • Drought monitoring
  • Flood assessment
  • Ocean mass change (sea level)
  • Terrestrial water storage
  • Earthquake mass redistribution

2. The Conceptual Model

Gravitational Acceleration

Newton’s law:

\[g = \frac{GM}{r^2}\]

Where:

  • $g$ = gravitational acceleration (m/s²)
  • $G$ = gravitational constant (6.674 × 10⁻¹¹ m³/kg/s²)
  • $M$ = mass (kg)
  • $r$ = distance (m)

Earth surface: $g \approx 9.81$ m/s²

Variations:

Latitude (centrifugal force): ±0.03 m/s²
Altitude (1 km): -0.003 m/s²
Mass anomalies: ±10⁻⁶ m/s² (1 μGal)

GRACE detects: ~10⁻⁸ m/s² (0.01 μGal)

Satellite Perturbations

Gravity anomaly → orbital velocity change

Two satellites in tandem:

Leading satellite over mass anomaly:

  • Accelerates (stronger pull)
  • Separation increases

Trailing satellite reaches anomaly:

  • Accelerates (catches up)
  • Separation decreases

Range rate:

\[\dot{\rho} = \frac{d\rho}{dt}\]

Where $\rho$ = inter-satellite distance (~220 km)

Measured: Range rate via K-band ranging (24 GHz)

Precision: 0.1 μm/s

Mass Anomaly Inversion

Gravitational potential:

\[V = \frac{GM}{r} + \text{anomalies}\]

Spherical harmonic expansion:

\[V = \frac{GM}{r} \sum_{l=0}^{\infty} \sum_{m=0}^{l} \left(\frac{a}{r}\right)^l P_{lm}(\sin\phi) (C_{lm}\cos m\lambda + S_{lm}\sin m\lambda)\]

Where:

  • $l$ = degree (spatial scale)
  • $m$ = order
  • $P_{lm}$ = associated Legendre polynomials
  • $C_{lm}$, $S_{lm}$ = Stokes coefficients
  • $a$ = Earth radius
  • $\phi$ = latitude
  • $\lambda$ = longitude

Truncation: $l_{\max} \approx 60$ (GRACE resolution limit)

Degree 60: ~300 km wavelength

Water Equivalent Thickness

Mass anomaly to water depth:

\[\Delta h = \frac{\Delta \sigma}{\rho_w}\]

Where:

  • $\Delta h$ = equivalent water thickness (m)
  • $\Delta \sigma$ = surface density anomaly (kg/m²)
  • $\rho_w$ = 1000 kg/m³

From gravity:

\[\Delta \sigma = \frac{a}{3} \sum_{l,m} \frac{2l+1}{1+k_l} \Delta C_{lm} Y_{lm}\]

Where $k_l$ = load Love number (elastic deformation).

3. Building the Mathematical Model

Satellite Acceleration

Gravity gradient along track:

\[\frac{\partial g}{\partial x} = -\frac{2GM}{r^3} + \text{anomaly gradient}\]

Differential acceleration:

\[\Delta a = \frac{\partial g}{\partial x} \times \Delta x\]

Where $\Delta x$ = satellite separation (~220 km)

Range rate change:

\[\ddot{\rho} = \Delta a\]

Integrate:

\[\Delta \rho(t) = \int_0^t \int_0^{t'} \Delta a \, dt' \, dt\]

Observed: Range vs predicted (from baseline gravity model)

Residual: Indicates mass change

Degree Variance

Power at each degree:

\[\sigma_l^2 = \sum_{m=0}^{l} (C_{lm}^2 + S_{lm}^2)\]

Kaula’s rule:

Expected variance:

\[\sigma_l \propto l^{-2}\]

GRACE measurement error:

Increases rapidly with degree:

\[\epsilon_l \propto l^2\]

Filtering required:

Low-pass (smooth) to reduce noise.

Gaussian filter:

\[W_l = e^{-l(l+1) b^2 / 2}\]

Where $b$ = smoothing radius (typically 300-500 km)

Trade-off: Noise reduction vs spatial resolution

Trend Estimation

Time series at location:

\[\Delta h(t) = a + b \times t + \sum A_i \cos(\omega_i t + \phi_i) + \varepsilon\]

Where:

  • $a$ = offset
  • $b$ = linear trend (mass change rate)
  • $A_i$ = seasonal amplitudes
  • $\omega_i$ = annual, semi-annual frequencies
  • $\varepsilon$ = noise

Least squares fit:

Solve for parameters.

Uncertainty:

Accounts for temporal correlation in residuals.

4. Worked Example by Hand

Problem: Calculate water storage change from GRACE.

Observations:

Region: California Central Valley (120°W, 37°N, radius 200 km)

GRACE data (simplified):

Month 1 (Jan 2022): Gravity anomaly = +50 μGal
Month 13 (Jan 2023): Gravity anomaly = -30 μGal

Change: -80 μGal

Calculate equivalent water thickness change.

Solution

Step 1: Convert gravity to mass

\[\Delta g = 2\pi G \Delta \sigma\]

Where $\Delta \sigma$ = surface density change (kg/m²)

\[\Delta \sigma = \frac{\Delta g}{2\pi G}\]

Units:

1 μGal = 10⁻⁸ m/s²

\[\Delta \sigma = \frac{-80 \times 10^{-8}}{2\pi \times 6.674 \times 10^{-11}}\] \[= \frac{-8 \times 10^{-7}}{4.19 \times 10^{-10}} = -1910 \text{ kg/m}^2\]

Step 2: Convert to water depth

\[\Delta h = \frac{-1910}{1000} = -1.91 \text{ m}\]

1.9 meters of water loss!

Step 3: Total volume

Area of region (circle, r = 200 km):

\[A = \pi r^2 = \pi \times (200000)^2 = 1.26 \times 10^{11} \text{ m}^2\]

Volume change:

\[V = A \times \Delta h = 1.26 \times 10^{11} \times (-1.91) = -2.4 \times 10^{11} \text{ m}^3\]

= 240 km³ loss

Step 4: Interpretation

240 cubic kilometers of water lost in one year.

Causes:

  • Groundwater extraction
  • Below-average precipitation
  • Snow deficit
  • Soil moisture depletion

This is severe drought (typical seasonal variation: ±50 km³)

5. Computational Implementation

Below is an interactive GRACE data simulator.

Linear trend: -- cm/year

Seasonal amplitude: -- cm

Total change: -- m

Mass change rate: -- Gt/year

Observations:

  • Greenland shows strong negative trend (ice loss: -280 Gt/year)
  • Amazon shows large seasonal variations (wet/dry seasons)
  • California shows drought signal with seasonal variation
  • Ganges shows groundwater depletion trend
  • Blue line: actual GRACE signal (trend + seasonal + noise)
  • Red dashed: linear trend component
  • Seasonal variations evident as annual oscillations

Key findings:

  • GRACE detects both long-term trends and seasonal cycles
  • Ice sheet mass loss clearly visible as negative trend
  • Groundwater depletion measurable in major aquifers
  • Seasonal water storage changes reach 10-25 cm equivalent

6. Interpretation

Ice Sheet Mass Balance

Greenland (2002-2023):

GRACE trend: -280 Gt/year average

Acceleration: -20 Gt/year² (increasing loss)

Contributions to sea level:

280 Gt/year ÷ (ocean area × density) = 0.8 mm/year

Regional patterns:

  • Southeast: Largest losses
  • Northwest: Moderate losses
  • Interior: Slight gains (snow accumulation)

Antarctica:

GRACE trend: -150 Gt/year

Variations:

  • West Antarctica: -160 Gt/year (marine ice sheet collapse)
  • East Antarctica: +10 Gt/year (slight snow increase)
  • Antarctic Peninsula: -20 Gt/year

Combined: ~1 mm/year sea level rise from ice sheets

Groundwater Depletion

North India (Ganges-Brahmaputra):

GRACE: -4 cm/year water storage loss

Cause: Irrigation extraction > recharge

Volume: -54 km³/year

Unsustainable: Fossil aquifer depletion

California Central Valley:

2011-2015 drought:

  • GRACE: -15 cm/year peak loss
  • Groundwater contributed 60% of deficit
  • 50+ km³ cumulative loss

Recovery: 2017-2019 wet years partially replenished

Drought Monitoring

2010-2011 Amazon drought:

GRACE detected:

  • -15 cm water storage anomaly
  • Preceded vegetation stress (optical NDVI)
  • Early warning capability

Operational use:

  • USDA drought monitor integrates GRACE
  • NASA FLDAS (Famine Early Warning System)
  • Water resource planning

7. What Could Go Wrong?

Spatial Leakage

Smoothing spreads signal:

300 km Gaussian filter → adjacent regions contaminate

Example:

Greenland ice loss “leaks” to ocean, land nearby.

Correction:

Forward modelling:

  • Assume spatial pattern (coast concentration)
  • Apply GRACE processing
  • Compute gain factors
  • Amplify observed signal

Typical: 10-30% underestimate without correction

Glacial Isostatic Adjustment (GIA)

Ice age deglaciation:

Mantle still rebounding.

GIA vertical motion:

Up to 10 mm/year (Scandinavia, Canada)

GRACE sees mass change:

Cannot distinguish GIA from contemporary changes.

Correction:

GIA models (ICE-6G, etc.) based on ice history.

Subtract from GRACE signal.

Uncertainty: ±20-30% in some regions

Geocenter Motion

Earth’s center of mass moves:

Relative to crust surface.

Degree 1 (l=1) coefficients:

Not measured by GRACE (both satellites affected equally).

Estimate from:

  • Ocean models
  • Station position networks
  • Combination solutions

Impact: Small global (few mm), important for sea level budget

Earthquake Signals

Large earthquakes:

Redistribute mass (coseismic + postseismic).

2004 Sumatra M9.1:

GRACE detected:

  • Coseismic: -5 cm water equivalent (localized)
  • Postseismic: Years of relaxation

Challenge:

Separate earthquake from hydrologic signals.

Solution:

Model earthquake, subtract before hydrology analysis.

8. Extension: GRACE Follow-On

GRACE-FO (launched 2018):

Improvements:

  • Laser ranging (addition to microwave)
  • Better accelerometers
  • Continuous from GRACE

Laser ranging:

10-100× better precision than microwave.

But: Atmospheric scattering limits (clouds block laser).

Combined system:

Microwave for continuous tracking, laser for highest precision.

Future missions:

Mass Change mission (2028+):

  • Lower orbit (improve resolution)
  • Better instruments
  • Goal: 150 km resolution

Applications expansion:

  • Smaller aquifers
  • Individual drainage basins
  • Urban water use
  • Irrigation monitoring

9. Math Refresher: Spherical Harmonics

Basis Functions

On sphere:

\[Y_{lm}(\theta, \phi) = P_{lm}(\cos\theta) e^{im\phi}\]

Where:

  • $\theta$ = colatitude
  • $\phi$ = longitude
  • $P_{lm}$ = associated Legendre polynomial

Properties:

Orthogonal:

\[\int Y_{lm} Y_{l'm'}^* \, d\Omega = \delta_{ll'} \delta_{mm'}\]

Complete: Any function on sphere can be expanded.

Wavelength

Degree $l$ corresponds to wavelength:

\[\lambda = \frac{2\pi a}{l}\]

Where $a$ = Earth radius (6371 km)

Examples:

  • $l = 2$: ~20,000 km (Earth’s flattening)
  • $l = 10$: ~4,000 km (continents)
  • $l = 60$: ~670 km (GRACE resolution)
  • $l = 360$: ~110 km (EIGEN-6C4 model)

Summary

  • GRACE satellites measure gravity variations from mass redistribution via inter-satellite ranging
  • Gravitational acceleration changes of 10⁻⁸ m/s² detectable enabling water storage monitoring
  • Mass anomalies inverted from spherical harmonic coefficients of gravitational potential
  • Water equivalent thickness derived from surface density changes via gravity relationship
  • Temporal resolution monthly with spatial resolution approximately 300-400 km
  • Applications span ice sheet mass balance, groundwater depletion, drought monitoring
  • Greenland losing 280 Gt/year, Antarctica 150 Gt/year contributing to sea level rise
  • Challenges include spatial leakage, GIA correction, geocenter motion estimation
  • GRACE Follow-On continues measurements with improved laser ranging capability
  • Critical tool for water resources assessment and climate change monitoring

References