Is a glacial lake above my valley growing toward an outburst flood?

A glacial lake is a pond of meltwater that forms where a glacier has retreated. It is often held back by a loose natural wall of rock and rubble called a moraine, not solid bedrock. When that wall fails, the whole lake can drain at once into a wall of water rushing down the valley: a glacial lake outburst flood (GLOF). The most useful early-warning sign is simple. A lake that keeps getting bigger, especially faster and faster, is putting more pressure on the wall holding it in.

Which data to use, and how

Use this dataset: HLS Sentinel-2 (short name HLSS30). It gives you cloud-screened satellite pictures at 30 m per pixel from 2015 onward, all through one free NASA login. “HLS” means Harmonized Landsat Sentinel-2: NASA has already stitched two satellites’ images onto the same grid so you don’t have to. That makes it the easiest single-login place to start. Want sharper detail or a longer history? Sentinel-2 MSI L2A is finer (10 m, 2015+), and Landsat 8/9 C2 L2 reaches back to 1984. Add those later if you need them.

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

1. Pick your valley. Draw a small latitude/longitude box around the lake and choose the years you want to scan.
2. Grab one clear picture per year. Aim for the post-monsoon window (roughly October), when the sky is clearest and the seasonal snow that fools the detector is at its lowest.
3. Find the water. Compute a water index, a simple formula that compares how bright each pixel is in two colour bands so water (dark in near-infrared) stands out from rock, ice, and snow. The common one is NDWI = (green − near-infrared) / (green + near-infrared).
4. Outline the lake and add up its size. Keep only pixels above a cutoff (the threshold), then multiply the number of water pixels by the area each pixel covers (30 m × 30 m) to get the lake area in km².
5. Plot area vs year and watch the slope. Rising area is the warning sign. Rising and speeding up is the strongest one. You can also overlay terrain (the DEM below) to see how steep the ground is between the lake and the nearest village.

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

What you can find out

• Whether a glacial lake is growing year over year for any high-mountain box you choose.
• Whether that growth is speeding up. This is the strongest early-warning signal for an outburst flood.
• Which lake in a basin is expanding fastest when there are several, as a quick screening pass.
• The shape of the natural dam (the moraine and lake-surface elevation), using the Copernicus GLO-30 elevation model and ICESat-2 along-track height measurements.
• A rough sense of who is downstream and exposed: how steep and how far the drop is between the lake and the nearest settlement.

What it can’t tell you

• Whether a flood is about to happen. Growing area raises the risk but does not trigger the flood. The actual breach can come from a rock or ice avalanche splashing into the lake, water seeping through the wall, or the loose interior of the moraine giving way. None of that shows up in a picture of the lake’s surface.
• Lakes you can’t see from above. Ponds hidden under rock debris or beneath the glacier’s skin (debris-covered and supraglacial lakes) are easily missed, because the water index only catches open, bright water.
• What’s happening during the monsoon. From roughly June to September, thick cloud hides the ground for months. Seeing through cloud needs radar (e.g. Sentinel-1 SAR), which bounces signals off the surface regardless of weather.
• An official hazard rating. Saying how a dam would actually break and how far the water would travel needs field surveys and physics-based flood models. This method only flags candidates for an expert to look at.

Gotchas to watch for

• The water index gets fooled. Cloud shadow, mountain shadow, and fresh snow can all look like water, while a muddy (turbid) or partly frozen lake can slip through unseen. Any false patch quietly inflates or shrinks your area number, so look at the outlines on the actual image every year (visual QA) before trusting them.
• The cutoff is not a universal number. Copying someone else’s threshold mislabels pixels, because the value that separates water from not-water depends on the local lake colour and snow. Start somewhere near 0 to 0.2 for NDWI, then tune it per scene against what you can see in the image.
• The monsoon leaves a gap. From June to September thick cloud blanks out the optical view for months, so part of the year may have no clear picture at all. Lean on the Landsat long record to fill in history, and add Sentinel-1 radar when you need all-weather monitoring.
• This is screening, not a verdict. A growth flag means “an expert should look at this lake,” not “a flood is coming.” Treat it as a shortlist, never as a hazard assessment.
• Different satellites must line up. Comparing area across sensors only works when they sit on the same map grid (consistent reprojection / co-registration); even a small misalignment makes a lake look like it changed when it didn’t. HLS already harmonizes Landsat and Sentinel-2 for you, which is another reason to start there.

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 HLS, open the green/NIR bands, compute NDWI, threshold to a lake mask, total the wetted area per year, fit the growth trend, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, rasterio
from rasterio.warp import transform_bounds
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

aoi = (78.5, 30.3, 80.5, 31.3)   # your valley as (W, S, E, N) — here, Uttarakhand Himalaya
years = range(2016, 2026)

# 1. One lake-area number per year. For each year grab HLS scenes in the
#    post-monsoon window (least cloud, low fresh snow), pick the least-cloudy one,
#    read green (B03) + NIR (B05), compute NDWI, threshold, and sum wetted pixels.
def lake_area_km2(yr):
    grans = earthaccess.search_data(
        short_name="HLSS30", bounding_box=aoi,
        temporal=(f"{yr}-10-01", f"{yr}-11-15"))
    if not grans:
        return np.nan
    # least-cloudy granule wins (HLS carries a scene cloud-cover attribute)
    def cloud(g):
        try:
            return float(g["umm"]["AdditionalAttributes"][0]  # CLOUD_COVERAGE
                         and next(a["Values"][0] for a in g["umm"]["AdditionalAttributes"]
                                  if a["Name"] == "CLOUD_COVERAGE"))
        except Exception:
            return 100.0
    best = min(grans, key=cloud)
    links = [u for u in best.data_links() if u.endswith((".B03.tif", ".B05.tif"))]
    files = earthaccess.open(sorted(links))  # B03 then B05
    with rasterio.open(files[0]) as g, rasterio.open(files[1]) as n:
        # clip both bands to the AOI in the scene's UTM grid
        w, s, e, no = transform_bounds("EPSG:4326", g.crs, *aoi)
        win = g.window(w, s, e, no)
        green = g.read(1, window=win).astype("float32")
        nir = n.read(1, window=n.window(w, s, e, no)).astype("float32")
        px_m = abs(g.transform.a)                  # ~30 m
    ndwi = (green - nir) / (green + nir + 1e-6)
    lake_mask = ndwi > 0.2                          # threshold — tune & QA visually per scene
    return float(lake_mask.sum()) * px_m * px_m / 1e6

area = np.array([lake_area_km2(y) for y in years])
yrs = np.array(list(years))[~np.isnan(area)]
area = area[~np.isnan(area)]

# 2. Growth trend (slope) and acceleration (curvature) — the GLOF precursors.
slope = np.polyfit(yrs, area, 1)[0]            # km^2 / yr
accel = np.polyfit(yrs, area, 2)[0]            # curvature: + => speeding up
growing = slope > 0
print(f"Lake area {'GROWING' if growing else 'stable/shrinking'}: "
      f"{slope:+.4f} km2/yr, curvature {accel:+.5f} "
      f"({'ACCELERATING — flag for expert review' if growing and accel > 0 else 'no acceleration'})")

# 3. Plot area vs year with the fitted trend.
plt.plot(yrs, area, "o-", label="lake area")
plt.plot(yrs, np.polyval(np.polyfit(yrs, area, 2), yrs), "--", c="r", label="quadratic fit")
plt.ylabel("lake area (km$^2$)"); plt.xlabel("year"); plt.legend(); plt.tight_layout(); plt.show()

# Before you trust it: eyeball each year's mask on the actual scene and cross-check the
# flagged lake against a published inventory (e.g. ICIMOD) — NDWI confuses shadow,
# turbid water, and snow for lake water (see the gotchas above). Screening, not a verdict.

Where the data comes from

This question pulls from several archives, but you only need one free NASA login (Earthdata Login) to begin. The starter dataset, HLS, comes from NASA’s LP DAAC and already merges Landsat and Sentinel-2 onto one grid, which is why it gives the fastest single-login start. The deeper layers live elsewhere: full-resolution Sentinel-2 and the Copernicus GLO-30 elevation model from ESA/Copernicus, the Landsat long record from USGS, and ICESat-2 elevation from NASA’s NSIDC DAAC. Each extra source can mean a separate login, so begin with HLS and add the others only as you need them.

Sources

• HLS (Harmonized Landsat Sentinel-2): https://hls.gsfc.nasa.gov/
• Copernicus DEM GLO-30: https://spacedata.copernicus.eu/collections/copernicus-digital-elevation-model
• ICESat-2 ATL06: https://nsidc.org/data/atl06
• ICIMOD glacial lakes / GLOF resources: https://www.icimod.org/
• NDWI water index (McFeeters 1996); MNDWI (Xu 2006) — see method notes before fixing thresholds