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

The Arctic is warming faster than anywhere else on Earth, and the tundra 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. NASA’s MOD11A2 gives a 1 km, 8-day map of how hot the land surface actually gets — not the air temperature a weather station reports, but the skin temperature of the soil, moss, and shrubs that sits directly above the frozen ground. Watch it across summers and you can see the surface heat that drives permafrost thaw.

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 and dark tundra surfaces can run far hotter than the air above them, which is exactly the heat that reaches down toward the permafrost table.

What you can answer

• How hot the tundra surface gets each summer, every 8 days back to 2000 — 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
• 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 is stored
• Whether the land is greening — overlay MOD13A2 NDVI for the same pixels to track shrub expansion and tundra “greening” alongside the warming
• Active-layer wetness context — SMAP soil moisture shows how wet the seasonally thawed surface layer is, which shapes both heat flow and plant growth

What you can NOT answer with these datasets alone

• Whether permafrost is actually thawing — MOD11A2 sees only the surface skin temperature. Permafrost is metres down; confirming thaw needs borehole temperatures or active-layer probes.
• How deep the active layer goes — surface temperature drives thaw but does not measure thaw depth; that requires ground sensors (e.g. CALM sites) or InSAR ground-subsidence analysis
• Cloud-free coverage every period — MODIS LST is optical/thermal; persistent Arctic cloud and the long polar night leave gaps, so winter and overcast composites can be sparse or missing
• Greening cause — NDVI rising could be more shrubs, longer growing seasons, or snow-timing shifts; the index alone does not separate them
• Sub-kilometre patchiness — tundra is a mosaic of polygons, ponds, and ridges; 1 km blends a warm dry ridge and a cool wet trough into one pixel

Code template (Python, cloud-direct)

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
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
g = earthaccess.search_data(short_name="MOD11A2",
                            temporal=("2023-07-01", "2023-07-31"),
                            bounding_box=(W, S, E, N))

# a bbox search returns multiple sinusoidal tiles; keep the one over the AOI
target = [x for x in g if "h12v02" in str(x)]
hdf = str(earthaccess.download(target[:1], local_path="./data")[0])

raw = SD(hdf, SDC.READ).select("LST_Day_1km").get().astype("float64")  # 1200x1200

# reconstruct lon/lat for MODIS sinusoidal tile h12 v02 to mask the AOI
R, T = 6371007.181, 1111950.519667
xmin, ymax, h, v = -20015109.354, 10007554.677, 12, 2
nrow, ncol = raw.shape
ii, jj = np.meshgrid(np.arange(ncol), np.arange(nrow))
px = T / ncol
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))

aoi = (lon >= W) & (lon <= E) & (lat >= S) & (lat <= N)
raw[raw == 0] = np.nan                       # 0 is the fill value
lst_c = raw * 0.02 - 273.15                   # raw -> kelvin -> °C
lst_c[~aoi] = np.nan
print("mean daytime LST:", round(float(np.nanmean(lst_c)), 2), "°C")

# greening: repeat with short_name="MOD13A2" and read the NDVI layer (scale 1e-4)
# trend: loop summers 2000..present, average July composites, fit a line