What’s happening with groundwater in my region?

Groundwater is the water sitting underground in the gaps of soil and rock, the stuff wells pump up for farms and cities. You can’t see it from space directly, but you can weigh it. When a region gains or loses water, the ground gets very slightly heavier or lighter, and a pair of NASA satellites measures that change in gravity from orbit. Track that weight over years and you can watch a region’s water (including the hidden groundwater) slowly drain away or refill.

Which data to use, and how

Use this dataset: GRACE-FO mascons (short name TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4). It gives you total water storage, all the water in a region added up (groundwater + soil moisture + snow + surface water), as a monthly map from 2002 to today. “Total water storage anomaly” (TWS) just means how much above or below the long-term normal the water is, in millimetres. Start here because one dataset already contains the whole signal.

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. Make it big (at least ~300 km across): these satellites only see large areas, so a small box won’t work. Pick your years too.
2. One number per month. Average the total-water-storage value inside your box for each month, giving you a monthly time series of how the region’s water rises and falls.
3. Subtract everything that isn’t groundwater. Total water storage includes soil moisture (from SMAP root-zone), snow, and surface water. Take those away and what’s left is the groundwater storage anomaly (GWSA), the key product.
4. Fit a trend. Draw the long-term line through the groundwater series. Down = the region is losing groundwater (depletion), up = it’s refilling, flat = stable. Pair it with rainfall (GPM IMERG) to see if dry years explain the dips.

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

• Whether your region’s total water is rising or falling over 20+ years, from GRACE+GRACE-FO.
• Surface and root-zone soil moisture trends from SMAP (2015 onward), the water in the top metre of soil.
• The groundwater part on its own: total water storage minus soil moisture, snow, and surface water. This is the headline number.
• Drought severity: how far below the 20-year average a given year sits.
• Irrigation-driven depletion in big farming aquifers — famous examples include California’s Central Valley, the Indo-Gangetic plain, the North China Plain, and the Middle East / North Africa (MENA).

What it can’t tell you

• Your individual well’s water level. The satellites see only big areas (~300 km across, the size of a small country). They can’t resolve a single well or a single farm.
• How much water a specific aquifer holds. Turning the storage change into an actual volume for one named aquifer needs extra ground data (geology, well logs) that the satellite doesn’t have.
• What happens next. Forecasting future groundwater needs a hydrological (water-cycle) computer model, not just past measurements.
• Fine-grained contrasts inside the region. You can’t reliably compare two spots closer than ~300 km apart. Splitting the signal finer (“downscaling”) means adding model guesses on top.

Gotchas to watch for

• 300 km is a hard limit. GRACE genuinely cannot see anything smaller, and claiming a result about one well or one field is the most common way people overclaim with this data. Keep your region big and talk only about region-wide trends.
• The subtraction can go wrong in snowy places. To isolate groundwater you subtract the soil moisture and snow estimates, but those come from a model, and in snow-heavy regions the snow estimate is uncertain. That error leaks straight into your groundwater number, so what you call “groundwater” could secretly be snow error. Be cautious in snowy or mountainous regions, and check the snow assumption against a model like LIS or NLDAS.
• There’s a gap in the record (2017–2018). The first satellite (GRACE) died and its replacement (GRACE-FO) launched about a year later, leaving a ~1-year hole. Ignore it and you can draw a fake trend straight across the gap. Handle the two records separately and bridge the gap carefully rather than pretending it’s continuous.
• Two flavours of the data disagree slightly. GRACE comes in “mascon” and “spherical-harmonic” versions, two math recipes that give slightly different answers, and mixing them muddies a trend. Pick one (mascons are easier) and use it consistently.
• Satellites show the combined signal, not the cause. A drop could be drought (natural) or pumping (human), and the data alone can’t assign blame. Compare against rainfall and climate to separate “no rain fell” from “people pumped it out.”

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, open the gridded mascon record, area-average to a monthly total-water series, subtract SMAP soil moisture to isolate groundwater, fit the trend, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, xarray as xr
import pymannkendall as mk
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

# Region: California Central Valley — keep the box big (>=300 km) for GRACE
aoi = (-122, 35, -119, 39)            # (W, S, E, N)
window = ("2002-04-01", "2025-12-31")

def area_mean(da, name):
    """Cosine-latitude-weighted mean of a grid clipped to the AOI."""
    da = da.sel(lon=slice(aoi[0] % 360, aoi[2] % 360), lat=slice(aoi[1], aoi[3]))
    w = np.cos(np.deg2rad(da.lat))
    return da.weighted(w).mean(("lat", "lon")).rename(name)

# 1. GRACE+GRACE-FO mascons — one continuous gridded record (cm of water).
grace = earthaccess.search_data(
    short_name="TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.3_V4",
    bounding_box=aoi, temporal=window)
gds = xr.open_dataset(earthaccess.open(grace)[0])
gds = gds.assign_coords(lon=(gds.lon % 360)).sortby("lon")
tws = area_mean(gds["lwe_thickness"] * 10.0, "tws")   # cm -> mm of water
tws = tws.resample(time="1MS").mean()                  # monthly; gap stays as NaN

# 2. SMAP L4 root-zone soil moisture (0-1 m) -> mm water-equivalent.
smap = earthaccess.search_data(short_name="SPL4SMGP", bounding_box=aoi,
                               temporal=("2015-04-01", "2025-12-31"))
sm_list = []
for fh in earthaccess.open(smap):
    sm = xr.open_dataset(fh, group="Geophysical_Data")["sm_rootzone"]  # m3/m3
    sm_list.append(float(area_mean(sm, "sm").values) * 1000.0)         # x1m -> mm
sm_series = np.array(sm_list)
sm_anom = sm_series - np.nanmean(sm_series)            # anomaly to match TWS frame

# 3. Rainfall context (mm/month) to tell drought from pumping.
imerg = earthaccess.search_data(short_name="GPM_3IMERGM", bounding_box=aoi,
                                temporal=window)
rain = xr.open_mfdataset([earthaccess.open(g) for g in imerg[:24]],
                         combine="nested", concat_dim="time")["precipitation"]
rain = area_mean(rain, "rain") * 24 * 30               # mm/hr -> ~mm/month

# 4. Isolate groundwater: GWSA = TWS - soil-moisture anomaly (snow ~ small here).
gws = tws.copy()
overlap = tws.sel(time=slice("2015-04-01", None))
gws.loc[dict(time=overlap.time)] = overlap.values - sm_anom[:overlap.size]

# 5. Robust trend on the groundwater series (NaNs across the 2017-18 gap dropped).
g = gws.dropna("time").values
r = mk.original_test(g)
trend = "depleting" if r.slope < 0 else "refilling"
print(f"Groundwater is {trend}: {r.slope*12:+.1f} mm/yr, p={r.p:.3f}",
      "(significant)" if r.h else "(not distinguishable from noise)")

# 6. Plot TWS, derived groundwater, and the fitted trend line.
plt.plot(tws.time, tws, c="0.7", label="total water storage")
plt.plot(gws.time, gws, c="C0", label="groundwater (GWSA)")
plt.axvspan(np.datetime64("2017-07"), np.datetime64("2018-06"),
            color="r", alpha=0.1, label="GRACE gap")
plt.ylabel("anomaly (mm)"); plt.legend(); plt.tight_layout(); plt.show()

# Before you trust it: compare the trend against in-situ well levels (USGS/DWR)
# and check the snow assumption against LIS/NLDAS (see the gotchas above).

Where the data comes from

Everything comes from NASA through one free login (Earthdata Login), even though the pieces live in different NASA archives: the gravity/total-water data (GRACE-FO) comes from PO.DAAC, the soil moisture (SMAP) from the NSIDC archive, the rainfall (IMERG) from GES DISC, and the optional land-model output (LIS) from GHRC. Four archives, one password.

Sources

• GRACE Tellus + groundwater app: https://grace.jpl.nasa.gov/applications/groundwater/
• SMAP L4 user guide: https://nsidc.org/data/spl4smau
• NASA NLDAS: https://ldas.gsfc.nasa.gov/nldas/