Is the ground under my city sinking — and is groundwater pumping to blame?

When an aquifer is pumped faster than it refills, the clay layers between the sand compact and the land surface slowly drops — sometimes by feet over a few dry years. Two NASA-supported records let you see both halves of the story: OPERA DISP-S1 turns Sentinel-1 radar into maps of how much the ground itself has moved, and GRACE-FO weighs the water underground by tracking tiny changes in Earth’s gravity. Together they connect sinking ground to vanishing groundwater.

Verified locally. Over California’s San Joaquin Valley AOI (36–37 °N, 120.6–119.6 °W), the GRACE-FO mascon lwe_thickness field shows the stored water mass falling at about -1.34 cm/yr of liquid-water-equivalent across the full record (Apr 2002 – Apr 2026), with the 2023 mean sitting at -15.18 cm below the mascon baseline. That is a clear, multi-decade groundwater drawdown signal — the kind of mass loss that drives the famous subsidence in this valley. Note: direct, city-scale InSAR subsidence needs OPERA DISP-S1 / Sentinel-1 (we found 244 granules over this AOI for 2023, but they stream from ASF DAAC, not GRACE’s gravity grid).

What you can answer

• How fast stored water is disappearing under a region — GRACE/GRACE-FO lwe_thickness (liquid-water-equivalent in cm), the mascon gravity signal of total water mass change
• How much the ground has moved, at ~30 m — OPERA DISP-S1 displacement (line-of-sight, in metres), built from Sentinel-1 InSAR so you can map subsidence bowls across a city or farm district
• Whether sinking tracks the water — overlay the InSAR displacement time series on the GRACE groundwater trend; if they fall together, pumping is the prime suspect
• The seasonal vs. permanent part — distinguish elastic rebound (land that springs back when wells recharge in winter) from inelastic compaction (the part that never comes back)
• How many people live on the sinking ground — drape WorldPop population over the subsidence map to turn millimetres into exposed residents and infrastructure

What you can NOT answer with these datasets alone

• The cause, with certainty — subsidence can also come from tectonics, oil/gas extraction, or peat oxidation; GRACE + InSAR show correlation, not a signed confession from a single well
• Which aquifer layer is compacting — InSAR sees the surface; pinning compaction to a specific clay interval needs extensometers or well-bore data
• City-block detail from GRACE — the mascon grid is coarse (~3°, ~300 km); a small AOI samples essentially one cell, so GRACE gives the regional water trend, not your street
• Absolute vertical motion from one track — OPERA DISP-S1 is line-of-sight; separating true vertical sinking from horizontal motion needs ascending + descending passes
• Future collapse or well-failure risk — these records are observational history, not a geotechnical forecast of when infrastructure will fail

Code template (Python, cloud-direct)

Verified locally. The headline number below came from the GRACE-FO mascon fallback (TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4, variable lwe_thickness, units cm), because the OPERA DISP-S1 granules — though present (244 over this AOI in 2023) — stream from ASF DAAC and need the ASF download path rather than a plain cloud-direct read. Swap in OPERA for true InSAR subsidence once you have ASF access configured.

import os, re, earthaccess, xarray as xr, numpy as np, pandas as pd

# load Earthdata creds from .env without `source` (passwords can break the shell)
for line in open(".env"):
    m = re.match(r'\s*(?:export\s+)?([A-Z0-9_]+)\s*=\s*(.*)\s*$', line)
    if m: os.environ.setdefault(m.group(1), m.group(2).strip().strip('"').strip("'"))
earthaccess.login(strategy="environment")   # free Earthdata Login

W, S, E, N = -120.6, 36.0, -119.6, 37.0     # San Joaquin Valley AOI

# --- groundwater mass trend from GRACE-FO (what we verified) ---
res = earthaccess.search_data(short_name="TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4")
ds  = xr.open_dataset(earthaccess.open(res[:1])[0])
lwe = ds["lwe_thickness"]                    # cm liquid-water-equivalent
lon = ds["lon"]
sel = dict(lat=slice(S, N),
           lon=slice(W % 360, E % 360) if float(lon.max()) > 180 else slice(W, E))
ts  = lwe.sel(**sel).mean(dim=["lat", "lon"]).load()

t    = pd.to_datetime(ts["time"].values)
days = (t - t[0]).days.values.astype(float)
m    = np.isfinite(ts.values)
slope = np.polyfit(days[m], ts.values[m], 1)[0] * 365.25
print("groundwater trend:", round(float(slope), 2), "cm/yr LWE")   # ~ -1.34

# --- direct InSAR subsidence (needs ASF DAAC access) ---
opera = earthaccess.search_data(short_name="OPERA_L3_DISP-S1_V1",
                                temporal=("2023-01-01", "2023-12-31"),
                                bounding_box=(W, S, E, N))
print("OPERA DISP-S1 granules found:", len(opera))   # 244
# files = earthaccess.download(opera[:1], "./opera")  # ASF download, then:
# disp  = xr.open_dataset(files[0])["displacement"]   # line-of-sight metres