Is the glacier above my valley thinning, and how fast?

A glacier is a slow river of ice. When it loses more ice to melting and to chunks breaking off the front than it gains from new snow, its surface drops. The top of the ice sits lower year after year. That drop is the clearest sign of a shrinking glacier, and a NASA satellite can measure it from space to within a few centimetres.

Which data to use, and how

Use this dataset: ATL06 (short name ATL06). It comes from NASA’s ICESat-2 satellite, which fires a green laser at the ground and times how long the pulse takes to bounce back. That travel time tells you the height of the surface. (Laser ranging from a satellite is called lidar.) ATL06 packages those returns into ready-to-use land-ice surface heights along narrow strips called ground tracks, the lines the satellite traces as it flies over. ICESat-2 has flown the same lines roughly every 91 days since 2018, so you can revisit one strip of a glacier and check whether it sits lower than it did before.

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

1. Pick your glacier: a small latitude/longitude box around it, and a stretch of time. The example uses the Columbia Glacier region of Alaska (60.3–60.7 °N).
2. Get the surface heights. Download the ATL06 files that cross your box and pull out the height value (h_li, in metres) for every laser footprint inside it.
3. Throw out the junk. Some footprints have no real reading (a placeholder “fill value” above 3×10³⁸). Drop those, and keep only the points that fall inside your box.
4. Compare across years. Find the same ground track from an earlier year, line the two passes up footprint-by-footprint, and subtract. The height difference is how much the surface dropped. Divide by the years between passes to get a thinning rate in metres per year.

Verified locally. For the Columbia Glacier region of Alaska, the ATL06 height field returned 20,760 valid surface-height measurements across 8 ICESat-2 passes in the 2023 melt season (March–September), spanning elevations from sea level up to 1,702 m with a median of 191 m. Note that thinning itself is the difference between repeat passes — ATL06 gives you the heights; you subtract matched tracks across years to get the rate.

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

What you can find out

• The surface height along a glacier, to within centimetres, every time ICESat-2 flew over.
• Whether the surface dropped between two years. Line up a repeat ground track from, say, 2019 and 2024 and subtract the heights to get a thinning amount in metres.
• A thinning rate. Divide that height change by the time between passes (metres per year).
• Where the loss is concentrated. Thinning is usually fastest near the terminus (the front end where the ice meets land or sea) and on the lower tongue, and slower up high where snow piles up.
• A long, consistent record. ICESat-2 has flown since 2018 with repeat tracks roughly every 91 days, giving many chances to revisit the same line.

What it can’t tell you

• A full thinning map of the whole glacier. The laser only measures along narrow tracks, not wall-to-wall. Between tracks you have blank gaps. To fill in the rest of the surface, pair it with a digital elevation model (a height map of the terrain). In this file’s list, that’s Copernicus DEM.
• How much ice or water was actually lost. A drop in height is not the same as a loss of mass. Turning metres of surface drop into tonnes of ice (or litres of water) needs an assumption about how dense the ice and snow are, plus the glacier’s area.
• Why it’s thinning. Melting, ice breaking off the front, and less snowfall all look the same in a height drop. Telling them apart needs climate data and where-the-front-sits data, not height alone.
• Where the ice front actually is. For the terminus position and how it retreats, you need optical pictures like HLSL30 (ordinary satellite imagery) alongside this.

Gotchas to watch for

• Clouds block the laser. A green laser can’t shoot through cloud, so cloudy passes simply have no data over your glacier and your coverage ends up patchy and uneven. Gather many passes across a season so the clear days fill in the gaps.
• Rough ice confuses the beam. Crevasses and jagged ice scatter the laser, so many footprints come back flagged as bad or missing and you won’t get a clean line everywhere. The product already marks these, and any decent workflow drops them.
• Fake “no reading” values look like real heights. Footprints with no valid measurement are stored as a huge placeholder number (a fill value above 3×10³⁸). Keep only heights below that threshold (the code does this with h < 3e38); forget to, and your average height will be wildly wrong.
• You must compare the same track. Subtracting two different tracks compares two different places, not change over time, so thinning only means something if the earlier and later passes trace the same line on the ground. Match passes by their track identifier (RGT / cycle) and line them up along-track before you subtract.

The real-data code

The Run it cell above runs this method on synthetic data (no login). Below is the whole recipe against the real archive: search two years of ATL06, open the beam heights, match the same ground track, subtract to get the surface drop, reduce it to a thinning rate, and plot it. It needs a free Earthdata Login.

import os, re, warnings
warnings.filterwarnings("ignore")
import earthaccess, h5py, numpy as np
import matplotlib.pyplot as plt

# 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 = -148.5, 60.3, -148.0, 60.7     # your glacier (Columbia, Alaska)
beams = ["gt1l","gt1r","gt2l","gt2r","gt3l","gt3r"]

# 1. Open ATL06 granules for one melt season and pull every good footprint,
#    keyed by its ground track (RGT) so we can match the SAME line across years.
def season_points(year):
    res = earthaccess.search_data(short_name="ATL06",
                                  temporal=(f"{year}-06-01", f"{year}-09-30"),
                                  bounding_box=(W, S, E, N))
    tracks = {}
    for fh in earthaccess.open(res):
        with h5py.File(fh, "r") as hf:
            rgt = int(hf["orbit_info/rgt"][0])           # ground-track id
            for b in beams:
                try:
                    h   = hf[f"{b}/land_ice_segments/h_li"][:]
                    lat = hf[f"{b}/land_ice_segments/latitude"][:]
                    lon = hf[f"{b}/land_ice_segments/longitude"][:]
                except KeyError:
                    continue
                m = (h < 3e38) & (lat>=S)&(lat<=N)&(lon>=W)&(lon<=E)   # drop fill values
                if m.sum():
                    tracks.setdefault(rgt, {"lat": [], "h": []})
                    tracks[rgt]["lat"].append(lat[m]); tracks[rgt]["h"].append(h[m])
    return {k: {"lat": np.concatenate(v["lat"]), "h": np.concatenate(v["h"])}
            for k, v in tracks.items()}

early, late = season_points(2019), season_points(2024)
yrs = 2024 - 2019

# 2. For each track flown in BOTH years, interpolate the later pass onto the
#    earlier pass's latitudes, subtract footprint-by-footprint, and stack.
lat_all, drop_all = [], []
for rgt in set(early) & set(late):
    a, b = early[rgt], late[rgt]
    o = np.argsort(b["lat"])
    h_late = np.interp(a["lat"], b["lat"][o], b["h"][o])  # match along-track
    lat_all.append(a["lat"]); drop_all.append(h_late - a["h"])
lat_all = np.concatenate(lat_all); drop_all = np.concatenate(drop_all)

# 3. Surface drop -> thinning rate (m/yr). Negative = the ice surface fell.
rate = drop_all / yrs
median_rate = np.median(rate)
verdict = "thinning" if median_rate < 0 else "thickening"
print(f"{len(set(early) & set(late))} repeat tracks; {rate.size} matched footprints")
print(f"Columbia Glacier surface is {verdict}: {median_rate:+.2f} m/yr (median, 2019->2024)")

# 4. Plot thinning rate against elevation: loss is fastest low near the terminus.
plt.scatter(lat_all, rate, s=4, alpha=0.3)
plt.axhline(0, c="k", lw=0.8); plt.axhline(median_rate, c="r", ls="--", label="median")
plt.xlabel("latitude (deg N, low->terminus)"); plt.ylabel("thinning rate (m/yr)")
plt.legend(); plt.tight_layout(); plt.show()

# Before you trust it: a difference only means change if both passes are the SAME
# track (matched on RGT, as above) — otherwise you compared two different places.

Where the data comes from

ATL06 comes from NASA through one free login (Earthdata Login): the ICESat-2 land-ice heights are served from the NSIDC archive (the National Snow and Ice Data Center). If you go further, the optical terminus imagery (HLSL30) and the terrain height map (Copernicus DEM) come from separate archives but through the same NASA login.