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

A glacier thins when it loses more ice to melt and calving than it gains from snowfall — and the clearest fingerprint is that its surface drops year after year. NASA’s ICESat-2 satellite fires a green laser at the ground and times the return — space-based lidar — to measure surface height to within a few centimetres. The ATL06 product turns those returns into land-ice surface heights along narrow ground tracks, so you can watch the same strip of a glacier and see whether it sat lower this year than last.

Verified locally. For the Columbia Glacier region of Alaska (60.3–60.7 °N), the ATL06 h_li field returned 20,760 valid surface-height measurements across 8 ICESat-2 passes in the 2023 melt season (Mar–Sep), spanning elevations from sea level up to 1,702 m with a median of 191 m. Thinning itself is the difference between repeat passes over the same track — ATL06 gives you the heights; you subtract matched tracks across years to get the rate.

What you can answer

• The surface elevation along a glacier, to within centimetres, every time ICESat-2 flew over — ATL06 h_li (metres; fill values above 3×10³⁸ must be dropped)
• Whether the surface dropped between two years — match a repeat ground track from, say, 2019 and 2024 and difference the heights to get a thinning amount in metres
• A thinning rate — divide that elevation change by the time between passes (m per year)
• Where the loss is concentrated — thinning is usually fastest near the terminus and on the lower tongue, slower up in the accumulation zone
• 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 you can NOT answer with these datasets alone

• A full thinning map of the whole glacier — ICESat-2 measures only along narrow laser tracks, not wall-to-wall. Between tracks you have gaps; pair with a digital elevation model for the full surface.
• Ice volume or mass loss directly — height change is not mass; converting to water needs a firn/ice density assumption and the glacier’s area
• Why it is thinning — melt versus calving versus reduced snowfall needs climate and terminus data, not elevation alone
• Cloudy-day or heavily-crevassed returns — clouds block the laser and rough ice scatters it; many segments are flagged or missing, so coverage is uneven
• The terminus position — for where the ice front sits and how it retreats, use optical imagery like HLSL30 alongside

Code template (Python, cloud-direct)

Verified locally. ATL06, variable land_ice_segments/h_li, is an HDF5 file with six beam groups (gt1l…gt3r). Search with a bounding box, drop fill values (> 3×10³⁸), and mask to your AOI. Repeat tracks across years, differenced, give thinning over time.

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

# 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)
res = earthaccess.search_data(short_name="ATL06",
                              temporal=("2023-03-01", "2023-09-30"),
                              bounding_box=(W, S, E, N))
files = earthaccess.download(res, local_path="/tmp/atl06")

beams = ["gt1l","gt1r","gt2l","gt2r","gt3l","gt3r"]
heights = []
for f in files:
    with h5py.File(f, "r") as hf:
        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)
            if m.sum():
                heights.append(h[m])

arr = np.concatenate(heights)
print("valid h_li over AOI:", arr.size)
print("elevation range: %.1f to %.1f m" % (arr.min(), arr.max()))

# thinning = repeat the search for an EARLIER year, match the same ground track
# (RGT / cycle), interpolate to common along-track positions, and subtract:
#   thinning_m_per_yr = (h_later - h_earlier) / years_between