Coseismic landslides are a major driver of landscape evolution and geohazards in active mountain belts, reshaping topography and influencing drainage divide migration. The debris they generate initiates sediment cascades, producing transient changes in sediment flux and river incision. However, accurately quantifying coseismic landslide erosion and deposition remains challenging. Following the 2016 Mw 7.8 Kaiko & strns;ura earthquake, we developed a 3D method to volumetrically map landslide erosion and debris distribution across hillslopes using a 2 m resolution vertical difference model over 6875 km2. Our results show that large coseismic landslides preferentially eroded rock mass from upper slopes, with most debris retained on-slope. Only 14% of total landslide volume was deposited off slope, with >= 4% within channels. In contrast, traditional 2D connectivity approaches estimated that up to 36% of landslides were channel-connected, with 89% of the total landslide volume available for fluvial transport. This new volumetric mapping demonstrates that the transfer of landslide debris from hillslopes to river systems is primarily governed by post earthquake remobilisation processes rather than immediate delivery following the earthquake.

This study investigates flooding and erosion impacts and human responses in Aoraki Mount Cook and Westland Tai Poutini national parks in Aotearoa New Zealand. These fast-eroding landscapes provide important test cases and insights for considering the public access dimensions of climate change. Our objectives were to explore and characterise the often-overlooked role of public access as a ubiquitous concern for protected areas and other area-based conservation approaches that facilitate connections between people and nature alongside their protective functions. We employed a mixed-methods approach including volunteered geographic information (VGI) from a park user survey (n = 273) and detailed case studies of change on two iconic mountaineering routes based on geospatial analyses of digital elevation models spanning 1986–2022. VGI data identified 36 adversely affected locations while 21% of respondents also identified beneficial aspects of recent landscape changes. Geophysical changes could be perceived differently by different stakeholders, illustrating the potential for competing demands on management responses. Impacts of rainfall-triggered erosion events were explored in case studies of damaged access infrastructure (e.g., roads, tracks, bridges). Adaptive responses resulted from formal or informal (park user-led) actions including re-routing, rebuilding, or abandonment of pre-existing infrastructure. Three widely transferable dimensions of public access management are identified: providing access that supports the core functions of protected areas; evaluating the impacts of both physical changes and human responses to them; and managing tensions between stakeholder preferences. Improved attention to the role of access is essential for effective climate change adaptation in parks and reserves.

Three-dimensional landscape changes were investigated in the Kitchener Avalanche Path, Aoraki/Mount Cook National Park, New Zealand, after an extreme storm in July 2022. The Path features an earthen diversion berm constructed in 2018 to mitigate the risk of avalanches to the adjacent Aoraki/Mount Cook Village. The berm performed as designed during a significant snow avalanche and related record-breaking winter rain event in 2022. However, damage to a section of the berm was observed from rain runoff during the storm. This study quantifies the impact of the storm on the Kitchener path runout zone, including the berm. Surveys conducted over five epochs (2008, 2018a, 2018b, 2022, and 2023) highlight the topographic changes from preconstruction through post-storm conditions. Using the derived digital elevation models, volume changes were estimated based on 2D cross-profiles and 3D DEM-differencing. Although significant erosion in the lower section of the berm was observed (-759 +/- 58 m3), minimal erosion along the inside face of the berm from avalanche impact was detected. Observed changes provide insights into the dynamic landscape, signalling that rain-induced erosion has significant impacts on engineered earthen protection structures. The future performance of the berm may require increasingly frequent repairs from extreme rain events causing further erosion.

Evaluating the influence of earthquakes on erosion, landscape evolution, and sediment-related hazards requires quantifying the number, location, and volume of co-seismic landslides following large earthquakes. Direct measurements of individual landslide volumes are difficult to accomplish and, as a result, rarely attempted. As an alternative, power-law scaling relationships between landslide volume and area—derived from global landslide inventories—are often used to estimate landslide volumes. The accuracy of the volume estimates derived from using these global scaling relationships remain untested and could be a source of significant uncertainty. Here, an ensemble method, adopting a pre- to post-event 2 by 2 m difference model covering an area of 6875 km2, has been developed to accurately estimate the volume of the > 30,000 mapped landslides triggered by the 2016 Mw 7.8 Kaikōura earthquake New Zealand. Our ensemble method estimates the total landslide volume for the Kaikōura earthquake to be 241 (+ 100/ − 69) M m3. This newly quantified total volume estimate of co-seismic landsliding was then compared to the total volumes estimated using globally and locally derived scaling relationships, which were found to underestimate the total landslide volume by up to 52

Background: Global canopy height models are becoming prolific yet require evaluation across New Zealand's diverse vegetation types to assess their accuracy and applicability. Accurate measurement of canopy height is crucial for estimating above-ground woody biomass, which is essential for modelling carbon emissions and sequestration in the context of climate change. These models generally rely on remote sensing data and machine learning techniques, with Light Detection and Ranging (LiDAR) technology commonly employed for precise measurement. Methods: This study validated the three latest global canopy height models, each provided at a different resolution: 30-metre, 10-metre, and 1-metre. We assessed the accuracy of the selected models by comparing them against canopy height estimates derived from local Airborne Laser Scanning (ALS) datasets, which served as our reference data. Eleven regions across New Zealand were selected based on ALS data availability, encompassing five vegetation and land cover types. Our methodology involved utilising and automating the processing of large New Zealand ALS datasets. To align resolutions for comparison, the reference canopy height was calculated by aggregating average or maximum heights at 10 and 30 m spatial resolution. Model performances were assessed using statistical metrics, including root-mean-square error (RMSE), bias, and R2. Results: Overall, all models exhibited relatively low R2 values, indicating limited capture of canopy height variability. The Potapov 30-metre model performed best with average aggregation in shorter vegetation. In contrast, the Lang 10-metre model showed improved accuracy with maximum aggregation, particularly in taller vegetation, but visual boundaries between different vegetation types were not as distinct. The Tolan 1-metre model provided a balanced approach, minimising biases in lower heights but underestimating taller canopies. Results highlight model-specific strengths for varying vegetation structures and the sensitivity of performances to aggregation methods applied to high-resolution reference ALS data. Conclusions: All three global canopy height models exhibit varied performance across New Zealand's vegetation types. The findings highlight the importance of vegetation-specific applications to optimise each global model's accuracy. Currently, these models are suitable for carbon accounting efforts as supplementary tools rather than replacements for existing methodologies.

Shrub encroachment into grassland ecosystems has been increasingly observed and documented worldwide in recent years. A grass–shrub transition can affect the diversity, abundance and functional integrity of grassland plant communities and understanding the drivers behind these processes is therefore crucial. While potential environmental drivers are often investigated, the role of spatial patterns of neighbouring shrub density in local shrub encroachment has been less well studied. The aim of this study is to investigate the relative role of neighbouring shrub density and topography as potential key drivers of shrub encroachment in a typical montane grassland ecosystem in New Zealand. We used the SPOT (Satellite Pour l’Observation) 6/7 multispectral imagery captured on one day in 2013 and in 2017 to calculate recent changes in shrub/grass cover during this period. Using the Normalised Difference Vegetation Index (NDVI), we classified the study area into grassland and shrubland and quantified the extent and change in these two land-cover types over the study period. We then investigated the relationships between changes in land cover and neighbourhood shrub density, elevation and aspect. Between 2013 and 2017, there was an overall shrubland increase of + 0.35

AbstractA distributed mass-balance model is used over a 10-year period for the re-analysis of a glaciological mass-balance time series obtained from Brewster Glacier, New Zealand. Mass-balance modelling reveals glaciological mass balance has been overestimated, with an average mass loss of −516 mm w.e. a−1 not captured by observations at the end of the ablation season, which represents 35% of the annual mass balance. While the average length of the accumulation season (199 days) remains longer than the ablation season (166 days), melting is not uncommon in the core part of the accumulation season, with 2–32% of total snowfall being melted. Refreezing of meltwater is also important, with 10% of surface and subsurface melt being refrozen in the present climate. Net radiation, driven primarily by net shortwave radiation, is the main contributor to melt energy, with melt variability mainly influenced by the turbulent heat fluxes, net longwave radiation and the heat flux from precipitation in the ablation season. Snowfalls in summer are an important moderator of melt, highlighting the critical role of the ice-albedo feedback and phase of precipitation on seasonal mass balance. A complete homogenisation of the long-term glaciological mass balance for Brewster Glacier is still required.

Annual mass balance is an important reflection of glacier status that is also very sensitive to climate fluctuations. However, there is no effective and universal albedo-based method for the reconstruction of annual mass balance due to the scarcity of field observations. Here, we present an improved albedo–mass balance (IAMB) method to estimate annual glacier surface mass balance series using remote sensing techniques. The averaged glacier-wide albedo derived with the MODImLab algorithm during the summer season provides an effective proxy of the annual mass change. Defined as the variation in the albedo as a function of elevation change, the altitude–albedo gradient (∂z/∂α) can be obtained from a glacier digital elevation model (DEM) and optical images. The Chhota Shigri glacier situated in the western Himalayas was selected to test and assess the accuracy of this method over the period from 2003 to 2014. Reconstructed annual mass budgets correlated well with those from the observed records, with an average difference and root mean square error (RMSE) of −0.75 mm w.e. a−1 and 274.91 mm w.e. a−1, respectively, indicating that the IAMB method holds promise for glacier mass change monitoring. This study provides a new technique for annual mass balance estimation that can be applied to glaciers with no or few mass balance observations.

This study proposes and evaluates the performance of two image matching algorithms applied to hillshades derived from consecutive high-resolution digital surface models (DSMs) to measure surface displacements on landslides. The method is applied on Te Horo, a slow-moving landslide in the Dart Valley of New Zealand’s Southern Alps/Kā Tiritiri o te Moana using pairs of hillshades derived from Airborne and Satellite Photogrammetric Mapping (APM/SPM) of imagery captured in 2018 and 2020. This novel approach uses the consistency of displacement predictions generated from multiple hillshade pairs to gauge match quality and mask unreliable displacement predictions. The study compares a widely used normalised cross-correlation algorithm (NCC) alongside an optical flow approach for image matching. The performances of both algorithms are assessed against manually derived displacements of prominent surface features. We demonstrate the effectiveness of the masking approach as well as the good performance of the optical flow algorithm in delivering dense, accurate displacement measurements efficiently, particularly when high resolution DSMs are available. The results show that the main translational body of the landslide was displaced at rates averaging 25–55 mm day−1 over the 2018–2020 period.

An exceptional July 2022 winter storm brought 550 mm of precipitation to the Southern Alps of New Zealand. A series of alpine mass movements occurred during the storm, including a widespread snow avalanche cycle, debris flows, and erosion from rain runoff. We detail the sequence of events in the Kitchener avalanche path. Here, two large snow avalanches were followed by a debris flow. Substantial erosion of deposition and the underlying alluvial fan were induced by runoff from over 300 mm of rain falling after the first avalanche. The Kitchener path saw the largest avalanche since 1986, testing the utility of a diversion berm constructed for a 1:100‐year event. Results from a unmanned aerial vehicle lidar survey and numerical modeling characterize the rain‐on‐snow hazard sequence. In particular, the rain‐on‐snow event occurred on a deep mid‐winter snowpack, offering insights into future hazards posed by increasingly frequent extreme alpine precipitation.

Evaluating the influence of earthquakes on erosion, landscape evolution and sediment-related hazards requires quantifying the volume and velocity of post-seismic sediment cascades. However, accurate estimates of post-earthquake sediment transfers remain rare. Following the 2016 MW7.8 Kaikōura earthquake in New Zealand, the volume of post-seismic erosion was quantified directly by measuring the ground surface change between 4 lidar surveys captured in 2016, 2017, 2019 and 2021 using the multiscale model-to-model cloud comparison (M3C2) algorithm. The lidar surveys covered the 62 km2 Hapuku and 66 km2 Kowhai river catchments within the Seaward Kaikōura Range, representing the two catchments with the highest density of co-seismic landsliding.The total co-seismic landslide source volume for the Hapuku Catchment was 30 ± 6 M m3,the catchment being dominated by a 17 M m3 rock avalanche which dammed the Hapuku River. In the 5 years after the earthquake a total of 10.60 ± 0.22 M m3 of sediment was post-seismically eroded (equivalent to ~26% of the co-seismic landslide debris volume when considering bulking of the landslide deposit). A total of 9.71 ± 0.23 M m3 of sediment was delivered to the riverbed resulting in considerable riverbed aggradation and 3.58 ± 0.28 M m3 was inferred to have been transported beyond the rangefront of the Seaward Kaikōura Range (equivalent to ~9% of the co-seismic landslide debris). The total co-seismic landslide source volume for the Kowhai Catchment was only 13 +4/-3 M m3. Over the 5 years 2.02 ± 0.10 M m3 of sediment was post-seismically eroded, equal to ~13% of the co-seismic landslide debris volume within the catchment. The volume delivered to the riverbed, 1.29 ± 0.10 M m3 and 0.85 ± 0.13 M m3 is presumed to have been transported beyond the rangefront (equivalent to ~5% of the co-seismic landslide debris).From these volumes, the rates at which the co-seismic landslide sediment was eroded from hillslopes, delivered off-slope to channels and exported from the range front were calculated. When projected, these rates of sediment conveyance suggest the volume of co-seismically generated sediment is likely to be evacuated from the rangefront within or close to the recurrence interval for ground motions equivalent to the Kaikōura earthquake. The Hapuku and Kowhai river catchments being examples of where co-seismic landsliding counterbalanced uplift.

Abstract An end of summer snowline (EOSS) photographic dataset for Aotearoa New Zealand contains over four decades of equilibrium line altitude (ELA) observations for more than 50 index glaciers. This dataset provides an opportunity to create a climatological ELA reference series that has several applications. Our work screened out EOSS sites that had low temporal coverage and also removed limited observations when the official survey did not take place. Snowline data from 41 of 50 glaciers in the EOSS dataset were retained and included in a normalised master snowline series that spans 1977–2020. Application of the regionally representative normalised master snowline series in monthly and seasonally resolved climate response function analyses showed consistently strong relationships with austral warm-season temperatures for land-based stations west of the Southern Alps and the central Tasman Sea. There is a trend towards higher regional snowlines since the 1990s that has been steepening in recent decades. If contemporary decadal normalised master snowline series trends are maintained, the average Southern Alps snowline elevation will be displaced at least 200 m higher than normal by the 2025–2034 decade. More frequent extremely high snowlines are expected to drive more extreme cumulative mass-balance losses that will reduce the glacierised area of Aotearoa New Zealand.

Grasslands in mountainous areas often show distinct responses in the timing of their growing season in relation to topographical variation. However, it is unclear which factors of topography (elevation, aspect, and slope) affect which phases of the timing in seasonal growth (start, peak, end, length of growing season). Here we investigated these relationships between topography and growing season timing in the three key grassland types in mountainous areas in South Island, New Zealand. From a near-daily NDVI (Normalized Difference Vegetation Index) dataset over a 16 year period (2001–2016), we extracted five annual land surface phenology indices: start, end, length, peak of the growing season and peak NDVI. Averages in these phenology indices were correlated with three topographical factors. The start of growing season occurred later by 7.1, 5.2 and 3.5 days per 100 m elevation in the three grassland types (Alpine, Tall Tussock and Low Producing grasslands). The end of the season occurred earlier by 1.8, 1.6 days and later by 0.4 days per 10-degree more south-facing (colder) aspects in the three grasslands. A longer growing season was observed at lower elevation and on north-facing (sunny) slopes in alpine grasslands. A later season peak occurred at higher elevation and on north-facing slopes in alpine grasslands and at higher elevations and on steeper slopes in non-alpine grasslands. Higher peak NDVI was detected at the lower elevation. Our results show that different facets of a landscape's topography affect different stages of a grassland's growing season, and these responses also differ between grassland types. This highlights the importance of considering all topographical features when relationships between the physical environment and biological responses are investigated.

A project in Aotearoa/New Zealand, is combining the use of high-quality DEMs from satellite photogrammetric mapping (SPM) with Lidar technologies to model hazards such as snow avalanches. The resulting topographic mapping can improve planning and preparedness to deal with mass movements of snow and debris in alpine regions.

Natural hazard models need accurate digital elevation models (DEMs) to simulate mass movements on real-world terrain. A variety of platforms (terrestrial, drones, aerial, satellite) and sensor technologies (photogrammetry, lidar, interferometric synthetic aperture radar) are used to generate DEMs at a range of spatial resolutions with varying accuracy. As the availability of high-resolution DEMs continues to increase and the cost to produce DEMs continues to fall, hazard modelers must often choose which DEM to use for their modeling. We use satellite photogrammetry and topographic lidar to generate high-resolution DEMs and test the sensitivity of the Rapid Mass Movement Simulation (RAMMS) software to the DEM source and spatial resolution when simulating a large and complex snow avalanche along Milford Road in Aotearoa/New Zealand. Holding the RAMMS parameters constant while adjusting the source and spatial resolution of the DEM reveals how differences in terrain representation between the satellite photogrammetry and topographic lidar DEMs (2 m spatial resolution) affect the reliability of the simulation estimates (e.g., maximum core velocity, powder pressure, runout length, final debris pattern). At the same time, coarser representations of the terrain (5 and 15 m spatial resolution) simulate avalanches that run too far and produce a powder cloud that is too large, though with lower maximum impact pressures, compared to the actual event. The complex nature of the alpine terrain in the avalanche path (steep, rough, rock faces, treeless) makes it a suitable location to specifically test the model sensitivity to digital surface models (DSMs) where both ground and above-ground features on the topography are included in the elevation model. Considering the nature of the snowpack in the path (warm, deep with a steep elevation gradient) lying on a bedrock surface and plunging over a cliff, RAMMS performed well in the challenging conditions when using the high-resolution 2 m lidar DSM, with 99 % of the simulated debris volume located in the documented debris area.

Context: It is important to understand the responses of alpine vegetation to recent anthropogenic climate change. The mountainous landscapes with high climatic heterogeneity are good locations to investigate the effects of microclimatic variation on alpine ecosystems. Objectives: a) To what degree do topographical factors (aspect and elevation) affect the timing of growing season in alpine grasslands? b) Are these topographical effects different on alpine and non-alpine grasslands? Methods: We extracted five annual growth phenology indices (Start, End, Length, Peak and Peak-NDVI) in alpine and non-alpine grasslands in the Clutha river catchment, New Zealand with a near-daily NDVI (Normalized Difference Vegetation Index) dataset through 16 years (2001-2016). The shifting rates of these phenology indices were quantified with two topographical factors (aspect and elevation). Results: The start of season was delayed by 7.5, 5.1 and 3.7 days per 100 m higher of elevation in three grassland types (Alpine, Tall Tussock and Low Producing) respectively, and the end of season was advanced by 1.7, 1.3 days and delayed by 0.3 days per 10-degree-south on slopes individually. The longer season length was observed at lower elevation and on north-facing (sunny) slopes. The later season peak occurred at higher elevation and on north-facing slopes. The lower peak NDVI was detected at the higher elevation. Conclusions: In the studied grasslands, aspect and elevation were correlated to different phenological indices, and they affect phenology independently. The topographical effects are more pronounced in alpine ecosystems at the elevation above 1,300 m than in non-alpine ecosystems at lower elevation.

© 2019 The Author(s). Published by IOP Publishing Ltd. During austral summer (DJF) 2017/18, the New Zealand region experienced an unprecedented coupled ocean-atmosphere heatwave, covering an area of 4 million km2. Regional average air temperature anomalies over land were +2.2 °C, and sea surface temperature anomalies reached +3.7 °C in the eastern Tasman Sea. This paper discusses the event, including atmospheric and oceanic drivers, the role of anthropogenic warming, and terrestrial and marine impacts. The heatwave was associated with very low wind speeds, reducing upper ocean mixing and allowing heat fluxes from the atmosphere to the ocean to cause substantial warming of the stratified surface layers of the Tasman Sea. The event persisted for the entire austral summer resulting in a 3.8 ± 0.6 km3 loss of glacier ice in the Southern Alps (the largest annual loss in records back to 1962), very early Sauvignon Blanc wine-grape maturation in Marlborough, and major species disruption in marine ecosystems. The dominant driver was positive Southern Annular Mode (SAM) conditions, with a smaller contribution from La Niña. The long-term trend towards positive SAM conditions, a result of stratospheric ozone depletion and greenhouse gas increase, is thought to have contributed through association with more frequent anticyclonic 'blocking' conditions in the New Zealand region and a more poleward average latitude for the Southern Ocean storm track. The unprecedented heatwave provides a good analogue for possible mean conditions in the late 21st century. The best match suggests this extreme summer may be typical of average New Zealand summer climate for 2081-2100, under the RCP4.5 or RCP6.0 scenario.

Seasonal snow dramatically alters surface-atmosphere exchanges of heat and moisture, particularly in the South Island of New Zealand. Despite this, detailed simulations of seasonal snow and comparisons to remotely sensed snow observations are lacking in New Zealand, partly due to uncertainties in near-surface meteorology in mountainous areas with sparse in-situ observations. Here we simulate seasonal snow cover and water storage across New Zealand on a 250 m x 250 m grid over a 3-year period (April 2017 to March 2020) with a simple snow model and near-surface meteorology extracted from the New Zealand Convective Scale Model (NZCSM). Simulations are validated against snow cover derived from MODerate Resolution Imaging Spectroradiometer (MODIS) satellite remote sensing observations. While NZCSM near-surface meteorology is generally colder and wetter than observations, we find that the spatial patterns and elevation distribution of simulated snow cover duration (SCD) agrees well with MODIS observations. Biases in SCD are found in areas, particularly the ranges east of the main divide, where simulated snow cover persists longer than observed. In contrast, alternate simulations using daily gridded observations of near-surface meteorology show a poor fit to MODIS SCD, with large areas having little or no simulated snow cover. The simulated seasonal cycle of New Zealand-wide total snow water storage shows a peak around 1 September equating to around half the long-term average monthly total rainfall for the South Island. The correspondence between simulated and MODIS snow covered area is found to be sensitive to the threshold used to define simulated snow cover, particular in early winter when widespread thin snow cover is common. To improve estimates of snow cover and water storage, future work should exploit new remote sensing products for validation and assimilation as well as disentangle uncertainties in snow model parameters and meteorological input using detailed meteorological and snow observations.

Snow depth has traditionally been estimated based on point measurements collected either manually or at automated weather stations. Point measurements, though, do not represent the high spatial variability in snow depths present in alpine terrain. Photogrammetric mapping techniques have progressed in recent years and are capable of accurately mapping snow depth in a spatially continuous manner, over larger areas and at various spatial resolutions. However, the strengths and weaknesses associated with specific platforms and photogrammetric techniques as well as the accuracy of the photogrammetric performance on snow surfaces have not yet been sufficiently investigated. Therefore, industry-standard photogrammetric platforms, including high-resolution satellite (Pléiades), airplane (Ultracam Eagle M3), unmanned aerial system (eBee+ RTK with SenseFly S.O.D.A. camera) and terrestrial (single lens reflex camera, Canon EOS 750D) platforms, were tested for snow depth mapping in the alpine Dischma valley (Switzerland) in spring 2018. Imagery was acquired with airborne and space-borne platforms over the entire valley, while unmanned aerial system (UAS) and terrestrial photogrammetric imagery was acquired over a subset of the valley. For independent validation of the photogrammetric products, snow depth was measured by probing as well as by using remote observations of fixed snow poles. When comparing snow depth maps with manual and snow pole measurements, the root mean square error (RMSE) values and the normalized median absolute deviation (NMAD) values were 0.52 and 0.47 m, respectively, for the satellite snow depth map, 0.17 and 0.17 m for the airplane snow depth map, and 0.16 and 0.11 m for the UAS snow depth map. The area covered by the terrestrial snow depth map only intersected with four manual measurements and did not generate statistically relevant measurements. When using the UAS snow depth map as a reference surface, the RMSE and NMAD values were 0.44 and 0.38 m for the satellite snow depth map, 0.12 and 0.11 m for the airplane snow depth map, and 0.21 and 0.19 m for the terrestrial snow depth map. When compared to the airplane dataset over a large part of the Dischma valley (40 km2), the snow depth map from the satellite yielded an RMSE value of 0.92 m and an NMAD value of 0.65 m. This study provides comparative measurements between photogrammetric platforms to evaluate their specific advantages and disadvantages for operational, spatially continuous snow depth mapping in alpine terrain over both small and large geographic areas.

We investigate the temporal dynamics of shifts in phenological responses of a range of key stages of the growing season in New Zealand’s three indigenous grassland types over the last 16 years (2001–2016). A near-daily Normalized Difference Vegetation Index (NDVI) time series from MODerate Resolution Imaging Spectroradiometer (MODIS) was used to extract five annual growth phenology indices, namely the Start, End, Length, Peak and Peak NDVI of a growing season. The start of the growing season advanced (i.e. happened earlier) by a median of 7.2, 6.0 and 8.8 days per decade in Alpine, Tall Tussock and Low Producing grassland, whereas the end of the season advanced by a median of 4.5, 0.4 and 0.4 days in the three types respectively. The length of growing season was extended by 3.2, 5.2 and 7.1 days per decade in these three grassland types. Over 86% of the investigated grassland areas showed an advancing (earlier) start of the growing season, and 74% of Alpine grassland showed a trend toward an earlier end of season. Over 63% of all grassland types showed an increase in growing season length. A trend toward earlier growing season peak and overall increasing NDVI in the three grassland types indicate a tendency for increasing vegetation vitality in grassland ecosystems in recent years. The start of growing season was correlated with atmospheric pressure (negatively) and precipitation (positively) changes in winter–spring months, while the timing of the season end is positively correlated with air temperature and solar radiation in summer–autumn months. Our study shows that different grassland types differ in magnitude – but not in direction – of their recent shifts in timing of key growing season stages with high-alpine grasslands showing the strongest response. This study highlights the usefulness of remote sensing for monitoring ecosystem-level phenological shifts over large areas and long time periods.