Is my river’s snowmelt flood coming earlier as the mountains warm?

A spring flood on a snow-fed river is really the mountains emptying out: the snowpack that built up all winter melts, runs downhill, and the river rises. As the climate warms, that melt-out tends to start sooner — so the question is whether the snow over your headwaters is disappearing earlier than it used to. NASA’s MOD10A1 gives a daily, 500 m map of snow cover from the MODIS sensor, so you can watch the white retreat off your basin every spring back to 2000.

Verified locally. For the Colorado River headwaters (39.0–39.8 °N), the MOD10A1 NDSI_Snow_Cover field read 42.1% of the area still snow-covered (NDSI 40–100) on 1 April 2024 — a single daily snapshot of how full the snowpack still was that morning. Track that fraction day by day across a spring and the date it crosses, say, 50% gives you a clean “melt-out timing” number you can compare year to year.

What you can answer

• How much of your headwaters is snow-covered on any given day, back to 2000 — the fraction of pixels with NDSI_Snow_Cover between 40 and 100
• When the snow melts out each spring — find the date the snow-covered fraction drops below a threshold (e.g. 50%), the “snow disappearance date”
• Whether melt-out is trending earlier — line up that date across 20+ springs and look for a drift toward earlier dates as the mountains warm
• How this spring compares to normal — overlay this year’s melt curve on the multi-year average
• Where rain is piling onto the melt — add GPM IMERG rainfall, since spring rain-on-snow can accelerate the flood
• Whether the river downstream actually rose — pair with SWOT reach river heights to connect the snowpack story to the channel

What you can NOT answer with these datasets alone

• How much water is in the snow — MOD10A1 sees snow area, not snow depth or snow water equivalent; a thin late cover and a deep one can look identical. For water volume you need SNOTEL, SNODAS, or airborne snow surveys.
• The actual flood peak or discharge — snow cover is the driver; the flow rate needs a stream gauge (e.g. USGS) — SWOT gives height, not a calibrated discharge.
• What happens under clouds — MODIS is optical, so cloudy days leave gaps; you must gap-fill or composite, and stormy melt periods are exactly when clouds are worst.
• Sub-pixel or forested snow — at 500 m a dense canopy hides snow underneath, and small patches blur; melt timing in thick forest is least reliable.
• Why it melted early — the data show that timing shifted, not whether warm air, dust on snow, or low winter accumulation caused it.

Code template (Python, cloud-direct)

Verified locally. MOD10A1, variable NDSI_Snow_Cover, is a daily HDF-EOS2/HDF4 tile in the MODIS sinusoidal grid. There is no global lat/lon index to slice, so download one granule, read the 2400×2400 array, and clip to your basin’s row/column window. Values 0–100 are percent snow; values above 100 are flags (cloud, water, fill) — keep only 0–100.

import os, re, math, tempfile, earthaccess, numpy as np
from pyhdf.SD import SD, SDC

# 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 = -106.9, 39.0, -105.9, 39.8     # your headwaters (Colorado River)
g = earthaccess.search_data(short_name="MOD10A1",
                            temporal=("2024-04-01", "2024-04-05"),
                            bounding_box=(W, S, E, N))
f = [str(p) for p in earthaccess.download(g[:1], tempfile.mkdtemp()) if str(p).endswith(".hdf")][0]
arr = SD(f, SDC.READ).select("NDSI_Snow_Cover").get()   # (2400, 2400) uint8

# clip the sinusoidal tile (here h09v05) to the AOI window
R, T, ncell = 6371007.181, 1111950.519667, 2400
cell = T / ncell
x0 = -20015109.354 + 9 * T          # tile h09 left edge (m)
y1 = 10007554.677 - 5 * T           # tile v05 top edge (m)
def rc(lon, lat):
    x = R * math.radians(lon) * math.cos(math.radians(lat))
    y = R * math.radians(lat)
    return (y1 - y) / cell, (x - x0) / cell
rows = [rc(lo, la)[0] for lo in (W, E) for la in (S, N)]
cols = [rc(lo, la)[1] for lo in (W, E) for la in (S, N)]
r0, r1 = max(0, int(min(rows))), min(ncell, int(max(rows)) + 1)
c0, c1 = max(0, int(min(cols))), min(ncell, int(max(cols)) + 1)
sub = arr[r0:r1, c0:c1]

snow  = (sub >= 40) & (sub <= 100)          # snow-covered pixels
valid = (sub >= 0)  & (sub <= 100)          # exclude cloud/fill flags
print("snow-covered fraction:", round(100 * snow.sum() / valid.sum(), 1), "%")

# melt-out timing: loop spring days, compute this fraction each day, and find the
# date it first drops below 50% — then compare that date across many springs