Is the Arctic tundra near me warming and greening as permafrost thaws?

The Arctic is warming faster than anywhere else on Earth. The tundra (the treeless, frozen-ground landscape near the poles) shows it two ways. The ground surface gets hotter in summer, and the thin skin of plants on top gets greener as shrubs creep north. What you actually measure here is land-surface temperature. That’s not the air temperature a weather station reports. It’s the skin temperature of the soil, moss, and shrubs sitting directly above the frozen ground. Bare, dark tundra can run far hotter than the air above it, and that surface heat is what reaches down toward the permafrost (ground that stays frozen year-round).

Which data to use, and how

Use this dataset: MOD11A2 (short name MOD11A2), NASA’s MODIS-Terra land-surface temperature product. It gives you a 1 km map every 8 days, back to 2000, long enough to compare summers across two decades. Start here because surface heat is what drives thaw, and the record is long and consistent.

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

1. Pick your tundra spot. Draw a small latitude/longitude box around it, and choose the summer months and years you care about.
2. Read the surface temperature. Pull the LST_Day_1km layer. It comes as raw counts, so scale it: raw × 0.02 gives kelvin, then subtract 273.15 to get °C.
3. Average over your box. Combine the valid 1 km pixels inside your area into one number per 8-day composite, then per summer.
4. Compare summers across years. Average the July–August composites year by year and fit a long-term line. Up means the surface is getting hotter; flat means no real change. You can also overlay greenness (MOD13A2 NDVI) on the same pixels to see if shrubs are spreading.

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 watch each step work.

Verified locally. For the North Slope tundra AOI (68–70 °N, around −149 °E), the MOD11A2 LST_Day_1km field for the 8-day composite starting 26 Jun 2023 averaged 21.37 °C daytime land-surface temperature across 20,565 valid 1 km pixels (raw counts × 0.02 → K, then − 273.15 → °C). That is the warm midsummer surface skin, not the air. Bare, dark tundra surfaces can run far hotter than the air above them, and that is the heat that reaches down toward the permafrost table.

What you can find out

• How hot the tundra surface gets each summer, every 8 days back to 2000, using MOD11A2 LST_Day_1km (raw × 0.02 → kelvin; subtract 273.15 for °C).
• Whether summers are trending warmer. Average the July–August composites year by year and fit a trend over two decades.
• The day-vs-night surface swing. Pair LST_Day_1km with LST_Night_1km to see how much the surface cools at night, a clue to how much heat the ground is storing.
• Whether the land is greening. Overlay MOD13A2 NDVI (a greenness index) on the same pixels to track shrub expansion and tundra greening alongside the warming.
• How wet the thawed surface layer is. SMAP soil moisture shows the wetness of the seasonally thawed top layer (the “active layer”), which shapes both heat flow and plant growth.

What it can’t tell you

• Whether permafrost is actually thawing. MOD11A2 sees only the surface skin temperature, but permafrost is metres down. Confirming thaw needs borehole temperatures or active-layer probes in the ground.
• How deep the active layer goes. Surface temperature drives thaw, but it doesn’t measure thaw depth. That requires ground sensors (e.g. CALM monitoring sites) or InSAR ground-subsidence analysis (radar that detects the ground sinking).
• Cloud-free coverage every period. MODIS land-surface temperature is optical/thermal, so it needs a clear view. Persistent Arctic cloud and the long polar night leave gaps, so winter and overcast composites can be sparse or missing.
• Why it’s greening. Rising NDVI could mean more shrubs, longer growing seasons, or shifts in snow timing. The index alone can’t separate those causes.
• Sub-kilometre detail. Tundra is a fine mosaic of polygons, ponds, and ridges. A single 1 km pixel blends a warm dry ridge and a cool wet trough into one number.

Gotchas to watch for

• Raw counts need scaling. The temperature layer ships as integers, not degrees, so you have to apply raw × 0.02 → kelvin, then − 273.15 → °C. Skip that and your “temperatures” land in the thousands. It’s the single most common beginner mistake with this product.
• Zero means “no data,” not “0 K.” The fill value for missing pixels is 0, so set those to NaN before you average and only let the valid pixels into the mean. Leave them in and they’ll drag your average toward absurd cold.
• The file comes in map tiles, and a box can straddle several. The data sits on the MODIS sinusoidal grid, cut into tiles, and a bounding-box search can return several at once. Pick the one covering your area (here h12v02 for the North Slope); grab the wrong tile and you’ll mask out your whole region by accident.
• Clouds and polar night punch holes. Because it’s a thermal sensor, cloudy or dark periods give thin or empty composites. Stick to summer composites and compare the same months across years (July vs July), where coverage is best.

The real-data code

The Run it cell above runs the method on synthetic data with no login. Below is the whole recipe against the real archive: search MOD11A2 tiles, open the HDF-EOS2 layer, mask the AOI, scale to °C, average each summer across two decades, fit a warming trend, overlay NDVI greening, and plot it. It needs a free Earthdata Login.

Verified locally. MOD11A2, variable LST_Day_1km, ships as HDF-EOS2/HDF4 tiles on the MODIS sinusoidal grid. A bounding-box search returns several tiles; pick the one covering your AOI (here h12v02 for the North Slope), read the layer, and scale raw counts × 0.02 → kelvin.

import os, re, earthaccess, numpy as np
import matplotlib.pyplot as plt
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 = -150, 68, -148, 70             # North Slope tundra, Alaska
TILE = "h12v02"                              # MODIS sinusoidal tile over the AOI
years = range(2003, 2024)

# Reconstruct lon/lat for sinusoidal tile h12 v02 (same grid for every 1200x1200 tile).
R, T = 6371007.181, 1111950.519667
xmin, ymax, h, v = -20015109.354, 10007554.677, 12, 2
ii, jj = np.meshgrid(np.arange(1200), np.arange(1200))
px = T / 1200
x = (xmin + h * T) + (ii + 0.5) * px
y = (ymax - v * T) - (jj + 0.5) * px
lat = np.degrees(y / R)
lon = np.degrees((x / R) / np.cos(y / R))
inbox = (lon >= W) & (lon <= E) & (lat >= S) & (lat <= N)

def aoi_mean(short_name, var, scale, offset=0.0, month_range=("07-01", "08-31")):
    """Download the AOI tile for one product/month-window and area-average to one number."""
    g = earthaccess.search_data(short_name=short_name, bounding_box=(W, S, E, N),
                                temporal=(f"{Y}-{month_range[0]}", f"{Y}-{month_range[1]}"))
    tiles = [x for x in g if TILE in str(x)]
    if not tiles: return np.nan
    hdf = str(earthaccess.download(tiles, local_path="./data")[0])
    raw = SD(hdf, SDC.READ).select(var).get().astype("float64")
    raw[raw == 0] = np.nan                   # 0 = fill value, not a real reading
    vals = raw * scale + offset
    vals[~inbox] = np.nan
    return float(np.nanmean(vals))           # area-average inside the AOI box

# Summer-mean July–August land-surface temperature, one number per year.
lst = []
for Y in years:
    lst.append(aoi_mean("MOD11A2", "LST_Day_1km", 0.02, -273.15))
lst = np.array(lst)
ok = ~np.isnan(lst)

# Fit a long-term warming line (°C per decade) over the summers we have.
yr = np.array(years)
slope, intercept = np.polyfit(yr[ok], lst[ok], 1)
trend = "warming" if slope > 0 else "cooling"
print(f"North Slope tundra summer LST is {trend}: {slope*10:+.2f} °C/decade "
      f"(mean {np.nanmean(lst):.1f} °C, {ok.sum()} summers)")

# Greening: same AOI, MOD13A2 NDVI (scale 1e-4) for the latest summer in the loop.
ndvi = aoi_mean("MOD13A2", "1 km 16 days NDVI", 1e-4)
print(f"latest-summer mean NDVI: {ndvi:.2f}  ({'greener' if ndvi > 0.3 else 'sparse'} tundra)")

# Plot the summers and the fitted trend.
plt.scatter(yr[ok], lst[ok], c="tab:red", label="summer LST")
plt.plot(yr[ok], intercept + slope * yr[ok], "k--", label=f"{slope*10:+.2f} °C/decade")
plt.ylabel("July–Aug daytime LST (°C)"); plt.xlabel("year")
plt.legend(); plt.tight_layout(); plt.show()

# Before you trust it: MOD11A2 is surface skin, not the permafrost metres below — confirm
# real thaw against borehole/active-layer probes (CALM sites), per the gotchas above.

Where the data comes from

Everything here comes from NASA through one free login (Earthdata Login). The surface temperature (MOD11A2) and greenness (MOD13A2) products both come from the LP DAAC archive, derived from the MODIS instrument on the Terra satellite. The optional wetness context, SMAP soil moisture, comes from a separate NASA satellite mission but uses the same single login.

Sources

• MOD11A2 land-surface temperature: https://lpdaac.usgs.gov/products/mod11a2v061/
• MOD13A2 vegetation indices (NDVI): https://lpdaac.usgs.gov/products/mod13a2v061/
• SMAP soil moisture: https://nsidc.org/data/spl3smp