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

Land subsidence is the slow sinking of the ground under your feet. One common cause: when an underground water store (an aquifer) is pumped faster than rain refills it, the clay layers squeezed between the sand collapse, and the whole land surface drops. In a few dry years that can be feet. Satellites watch both halves of this from orbit: how much the ground itself moves, and how much water is vanishing underground.

Which data to use, and how

Use this dataset: GRACE-FO mascon (short name TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4). GRACE-FO is a pair of satellites that “weighs” the water underground by sensing tiny changes in Earth’s gravity. When groundwater disappears, the ground gets slightly lighter, and the satellites feel it. The data come as a simple grid you can read straight from the cloud, which makes this the easiest starting point. (The other dataset here, OPERA DISP-S1, gives a much sharper ~30 m map of the actual sinking, but it streams from a different archive that needs extra download setup. Add it once you’re comfortable.)

Then do four steps (the runnable code further down does all of this — this is just the idea):

1. Pick your region. A latitude/longitude box around the city, aquifer, or farm district you care about, plus the years you want.
2. One water number per month. Read the GRACE-FO field called lwe_thickness (short for liquid-water-equivalent, in centimetres: how deep the water mass would be if you spread it out as a layer) and average it inside your box for each month.
3. Fit a trend. Draw the long-term line through that monthly series. A downward slope means stored water is steadily disappearing under the region.
4. Connect it to the sinking. Once you have OPERA DISP-S1 set up, overlay the ground-movement map on the water trend. If the ground falls as the water falls, pumping is the prime suspect.

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can see each step work.

What you can find out

• How fast stored water is disappearing under a region. The GRACE/GRACE-FO lwe_thickness signal, which is the mascon gravity measure of total water mass change.
• How much the ground has moved, at ~30 m. OPERA DISP-S1 displacement (how far the surface shifted along the radar’s line of sight, in metres), built from Sentinel-1 InSAR so you can map the “subsidence bowls” across a city or farm district.
• Whether the sinking tracks the water. Overlay the InSAR ground-movement series on the GRACE groundwater trend. If they fall together, pumping is the likely cause.
• The temporary vs. permanent part. Tell apart land that springs back when wells refill in winter (elastic rebound) from the sinking that never recovers (inelastic compaction).
• How many people live on the sinking ground. Drape WorldPop population data over the subsidence map to turn millimetres of sinking into a count of exposed residents and infrastructure.

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. (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 it can’t tell you

• The cause, for certain. Sinking can also come from tectonics (Earth’s crust shifting), oil/gas extraction, or peat soil drying out. GRACE + InSAR show that water loss and sinking happen together, but correlation isn’t a confession. Proving one specific well is to blame needs more data.
• Which underground layer is collapsing. InSAR only sees the surface. Pinning the squeeze to a particular clay layer needs in-ground instruments (extensometers) or well-bore measurements.
• City-block detail from GRACE. The gravity grid is coarse (~3°, roughly 300 km per cell), so a small box samples essentially one cell. GRACE gives the regional water trend, not your street.
• True up-and-down sinking from a single radar track. OPERA DISP-S1 measures motion along the satellite’s slanted line of sight, which mixes vertical sinking with sideways motion. Separating pure vertical needs both an ascending and a descending pass over the same spot.
• Future collapse or well-failure risk. These records are observational history (what already happened), not an engineering forecast of when infrastructure will fail.

Gotchas to watch for

• GRACE only sees big areas. Each gravity cell averages a ~300 km patch, and a tiny AOI just re-reads one fat pixel, so don’t expect neighbourhood detail. Treat the GRACE trend as “the whole region is losing water,” then turn to InSAR for the local picture.
• InSAR motion is slanted, not straight down. OPERA’s displacement is line-of-sight: motion measured along the radar beam’s diagonal, not straight vertical, so a raw number undercounts the true sinking. Combine an ascending and a descending track to recover the vertical, or just flag your value as a line-of-sight floor.
• OPERA lives in a different archive. The InSAR data stream from ASF DAAC (the Alaska Satellite Facility) and need that download path, not a plain cloud read like GRACE. That’s exactly why the headline number here came from the GRACE-FO fallback. Start with GRACE, configure ASF access, then swap OPERA in.
• Seasonal bounce isn’t permanent loss. Land often dips in summer (pumping season) and rebounds in winter (recharge), so comparing summer to winter can fake a “trend.” Compare the same season across years, or fit the full multi-year line.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search GRACE-FO, open the mascon grid, average lwe_thickness in your box, fit the groundwater trend, plot it, and confirm OPERA DISP-S1 InSAR exists for the sharper map. It needs a free Earthdata Login.

Verified locally. The headline number above came from the GRACE-FO mascon (TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4, variable lwe_thickness, units cm), because the OPERA DISP-S1 granules (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 earthaccess, xarray as xr, numpy as np, pandas as pd
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")          # free Earthdata Login

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

# 1. Search GRACE-FO mascons: one global gravity grid per month, 2002 -> present.
res = earthaccess.search_data(short_name="TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4")
print("GRACE-FO mascon granules:", len(res))

# 2. Open the mascon grid straight from the cloud and read liquid-water-equivalent.
ds  = xr.open_dataset(earthaccess.open(res[:1])[0])
lwe = ds["lwe_thickness"]                     # cm of water if spread as a flat layer
lon = ds["lon"]

# 3. Average inside the AOI for every month -> one water number per month.
#    The mascon grid runs 0-360 in longitude, so wrap the western box if needed.
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)

# 4. Fit a straight line through the monthly series -> slope in cm/yr.
#    A downward slope = stored water steadily vanishing under the region.
days  = (t - t[0]).days.values.astype(float)
good  = np.isfinite(ts.values)
slope = np.polyfit(days[good], ts.values[good], 1)[0] * 365.25
y2023 = float(ts.sel(time="2023").mean())     # 2023 mean vs the mascon baseline
verdict = "sinking water table" if slope < 0 else "recovering / stable"
print(f"groundwater is {verdict}: {slope:+.2f} cm/yr LWE; 2023 mean {y2023:+.2f} cm")  # ~ -1.34, -15.18

# 5. Plot the monthly LWE series with its fitted trend line.
plt.figure(figsize=(9, 3))
plt.plot(t, ts.values, lw=1, label="GRACE-FO LWE (cm)")
plt.plot(t[good], np.poly1d(np.polyfit(days[good], ts.values[good], 1))(days[good]),
         "r--", label=f"trend {slope:+.2f} cm/yr")
plt.axhline(0, c="k", lw=0.5); plt.ylabel("LWE anomaly (cm)")
plt.legend(); plt.title("San Joaquin Valley groundwater mass"); plt.tight_layout(); plt.show()

# 6. The sharper, direct check: OPERA DISP-S1 InSAR (~30 m line-of-sight sinking).
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, but they stream from ASF DAAC
# files = earthaccess.download(opera[:1], "./opera")  # configure ASF access, then:
# disp  = xr.open_dataset(files[0])["displacement"]   # line-of-sight metres

# Cross-check before trusting it: GRACE is coarse (~300 km/cell), so confirm the local
# sinking with the OPERA DISP-S1 InSAR map above — and recall line-of-sight undercounts
# true vertical until you combine an ascending and descending track (see the gotchas).

Where the data comes from

Both NASA records come through one free login (Earthdata Login), but from two different homes. The GRACE-FO gravity grid (your easy starting point) reads straight from the cloud. The OPERA DISP-S1 InSAR maps live at ASF DAAC (the Alaska Satellite Facility, NASA’s radar archive) and need its download path set up before you can pull them. The optional human-exposure layer (WorldPop population) is free and separate. So: two satellite missions, one login, one question. Is the ground sinking, and is the water why.