Subalpine treeline advance and shrubification monitoring under climate warming

Landsat's 50-year archive and Sentinel-2's 10 m resolution, combined with GEDI and ICESat-2 canopy-height profiles, let analysts track where alpine grassland ends and woody cover begins, and how that boundary has shifted across decades.

Sensors

Landsat 5/7/8/9 (TM / ETM+ / OLI): 30 m multispectral resolution, 16-day revisit per satellite, archive from 1984. The backbone for multi-decadal treeline trend detection; NDVI, NBR and tasselled-cap greenness derived from surface-reflectance products. Band 6 thermal on TM/ETM+ provides supplementary land-surface temperature context.
Sentinel-2 MSI (A and B combined): 10 m (visible/NIR) and 20 m (red-edge, SWIR) resolution, 5-day revisit at mid-latitudes with both satellites. Red-edge bands (B5, B6, B7 at 705–783 nm) improve discrimination of sparse shrub cover from bare soil and senescent grass, a distinction Landsat's broader bands struggle with.
GEDI (Global Ecosystem Dynamics Investigation): Spaceborne full-waveform lidar from the ISS, sampling a 25 m footprint at roughly 60 m cross-track spacing. Provides canopy-height estimates (RH95, RH100 metrics) that separate erect trees (typically >2 m) from prostrate or dwarf shrubs. Coverage limited to approximately 51.6° N/S latitude; many alpine zones fall within range.
ICESat-2 ATL08: Photon-counting lidar with a 17 m along-track footprint and 91-day repeat. ATL08 classifies ground and canopy photons to produce terrain-relative canopy height. Useful for validating GEDI-derived heights at the treeline ecotone and for polar or high-latitude sites outside GEDI's coverage.
Copernicus DEM (GLO-30): 30 m global digital elevation model derived from TanDEM-X radar. Not an imagery sensor, but an essential ancillary input: slope, aspect and hillshade layers computed from GLO-30 drive the topographic correction that removes shadow-induced underestimation of woody cover on pole-facing slopes.

What a 40-year NDVI trend actually tells you

The Landsat archive begins in 1984 with Thematic Mapper data. That is long enough to detect statistically significant greening trends at the pixel level across most alpine regions in the northern hemisphere. The standard approach stacks annual maximum NDVI composites, fits a Mann-Kendall trend test pixel by pixel, and maps where positive trends are concentrated along the upper forest belt. Sen's slope estimator then quantifies the rate of change in NDVI per year, giving a dimensioned, defensible figure rather than a qualitative impression.

NDVI alone cannot tell you whether a greening pixel contains a 4 m spruce or a 0.3 m crowberry mat. That distinction matters enormously for habitat assessment. The lidar products exist precisely to resolve it. Pairing a Landsat trend map with GEDI canopy-height samples lets you partition greening pixels into 'shrubification' and 'tree establishment' classes. Where GEDI footprints are sparse, Sentinel-2 red-edge indices (particularly the Chlorophyll Red-Edge index using B7 and B5) provide a continuous surface that correlates with shrub biomass and can be used to interpolate between lidar samples.

The shadow problem on pole-facing slopes

Topographic shadow is not a minor nuisance; on slopes steeper than 25° facing away from the sun, it can suppress apparent NDVI by 0.15 to 0.25 units relative to an equivalent sun-facing slope, according to published comparisons using the C-correction and Minnaert methods. Because treeline advance preferentially occurs on warmer, sun-facing aspects in most mountain systems, an uncorrected analysis will systematically underestimate woody cover on shadowed slopes and overestimate the asymmetry of advance. The result is a biased picture of how much alpine habitat has actually been lost.

The standard correction pipeline computes an illumination model from the DEM and solar geometry for each acquisition date, then applies either the C-correction or a semi-empirical Minnaert normalisation before deriving spectral indices. Neither method is perfect. The C-correction tends to overcorrect on very steep shadowed terrain, introducing artefacts. Minnaert requires a locally fitted k-parameter that varies by land-cover type. Both should be validated against field radiometry where possible. We flag this because any treeline product delivered without a documented topographic correction should be treated with caution, particularly if it claims to quantify pole-facing advance.

Tracking the closed-canopy boundary elevation

The treeline position is conventionally defined as the elevation at which closed canopy (typically >30% crown cover) persists for at least a decade. Mapping it from imagery requires a canopy-cover fraction product, not just a binary forest mask. Regression tree or random-forest models trained on high-resolution reference data and Sentinel-2 spectral bands can produce continuous fractional cover at 10 m, and the 30% contour can then be extracted per slope facet as an elevation value.

Repeat this extraction for each five-year epoch in the Landsat record and you have an elevation time series per aspect class. Published studies in the European Alps, Scandinavia and the Rocky Mountains have documented upslope advance rates ranging from roughly 1 to 4 metres of elevation per decade, though rates vary substantially by species, slope and precipitation regime. Those figures come from ground surveys and aerial photography cross-validated against satellite data; satellite-only estimates carry wider uncertainty, typically ±5 to 15 m in elevation depending on terrain complexity and image quality.

Shrub-cover fraction as the leading indicator

Shrubification often precedes tree establishment by decades. Dwarf willows, alders and ericaceous shrubs colonise fellfield and grassland at elevations where erect trees cannot yet survive, changing soil moisture, snow redistribution and invertebrate communities in ways that matter to conservation managers long before a canopy forms. Detecting this transition early requires the 10 m red-edge capability of Sentinel-2 rather than Landsat's 30 m bands.

A practical workflow derives a shrub-cover fraction map annually from Sentinel-2 summer composites, using spectral mixture analysis or a calibrated vegetation index. The fraction is then aggregated to 100 m grid cells for trend analysis, reducing noise from within-pixel heterogeneity. The honest limitation here is that prostrate shrubs under roughly 0.2 m height produce a spectral signature nearly indistinguishable from dense alpine grass in multispectral data. Field validation at the low end of the cover range is not optional; it is the only way to set a credible detection floor.

What the analysis cannot resolve, and why that matters

Several limits are worth stating plainly. Sentinel-2 revisit is 5 days in cloud-free conditions, but alpine terrain is frequently cloud-covered for weeks at a time. In practice, usable summer composites in many mountain ranges amount to 4 to 8 scenes per year, and in some years fewer. This compresses the phenological window available for shrub mapping and can introduce inter-annual artefacts if composite dates shift.

GEDI's sampling geometry means that in a 1 km² study area, you might have 20 to 80 footprints depending on latitude and acquisition history. That is sufficient for calibrating a continuous Sentinel-2 shrub-height product at landscape scale, but not for characterising individual patches smaller than roughly 0.5 ha. ICESat-2 adds density but not spatial continuity. Neither lidar system can detect woody stems below the canopy of taller vegetation, so understorey shrub encroachment beneath an existing forest canopy is invisible to these methods. Finally, attributing a detected trend to climate warming rather than changes in grazing pressure, fire history or land abandonment requires ancillary socioeconomic and management data that no satellite can supply.

Satellize's analytics pipeline for this use case runs on open Sentinel and Landsat archives, applies documented topographic correction, and delivers classified change layers with per-pixel uncertainty estimates. The approach is the same one underpinning the Tonga crop-estimation programme: open data, reproducible methods, auditable outputs.

Typical figures

Spatial resolution (optical)	10 m (Sentinel-2 visible/NIR), 20 m (Sentinel-2 red-edge/SWIR), 30 m (Landsat OLI/TM)
Revisit frequency	5 days (Sentinel-2A+B combined, cloud-free); 16 days per Landsat satellite; effective alpine summer composites typically 4–8 scenes per year
Canopy-height footprint	25 m (GEDI), 17 m along-track (ICESat-2 ATL08)
Archive depth	Landsat from 1984 (TM); Sentinel-2 from 2015; GEDI from 2019; ICESat-2 from 2018
Minimum detectable shrub-cover change	Approximately 5–10 percentage-point change in fractional cover at 100 m scale, subject to field validation; prostrate shrubs <0.2 m height are below reliable detection
Treeline elevation precision	±5–15 m in elevation, depending on terrain complexity and composite quality
Topographic correction	C-correction or Minnaert normalisation using Copernicus GLO-30 DEM (30 m); required on slopes >15° for reliable index derivation
Delivery formats	GeoTIFF change layers, GeoPackage vector contours, CSV elevation time series, PDF trend report
Latency (operational monitoring)	Annual update cycle typical; near-real-time flagging of anomalous greenness events possible within 10–15 days of cloud-free acquisition

Analytics Satellize can run

Multi-decadal NDVI trend map	Mann-Kendall trend test with Sen's slope on annual maximum NDVI composites from Landsat surface-reflectance archive (1984–present)	GeoTIFF raster of Sen's slope and significance mask; PDF summary report with elevation-stratified trend statistics
Treeline elevation time series	Fractional canopy-cover mapping from Sentinel-2 random-forest model; 30% cover contour extracted per aspect class per five-year epoch	CSV time series of median treeline elevation by aspect; GeoPackage vector contours per epoch
Shrub-cover fraction map (annual)	Spectral mixture analysis or calibrated red-edge vegetation index on Sentinel-2 summer composites; aggregated to 100 m grid	Annual GeoTIFF fractional-cover layers with per-pixel uncertainty band; change-detection layer comparing baseline and current period
Canopy-height classification (shrub vs. tree)	GEDI RH95 and ICESat-2 ATL08 canopy-height profiles used to train a height-class model applied to continuous Sentinel-2 spectral surface	GeoTIFF height-class map distinguishing prostrate shrub (<0.5 m), dwarf shrub (0.5–2 m) and erect tree (>2 m) zones
Topographic shadow correction report	Per-scene illumination modelling from GLO-30 DEM and solar geometry; C-correction applied; residual error quantified on north- and south-facing validation transects	Corrected surface-reflectance composites; QA report documenting correction method, validation transects and residual uncertainty by slope class
Aspect-stratified advance-rate summary	Elevation time series disaggregated by eight aspect classes; linear regression of treeline elevation against year per class	Summary table of advance rate (m elevation per decade) per aspect class with 95% confidence intervals; suitable for inclusion in conservation management plans

Sources

Who does the work

We can get this done for you. Satellize runs its own analyst desk and a strong science team. You do not buy a data feed and work out what it means; our people source the imagery, run the analysis described on this page, and hand you the answer with its confidence limits stated. Discuss this requirement.