Land cover changes across Greenland dominated by a doubling of vegetation in three decades

Land cover responses to climate change must be quantified for understanding Arctic climate, managing Arctic water resources, maintaining the health and livelihoods of Arctic societies and for sustainable economic development. This need is especially pressing in Greenland, where climate changes are amongst the most pronounced of anywhere in the Arctic. Ice loss from the Greenland Ice Sheet and from glaciers and ice caps has increased since the 1980s and consequently the proglacial parts of Greenland have expanded rapidly. Here we determine proglacial land cover changes at 30 m spatial resolution across Greenland during the last three decades. Besides the vastly decreased ice cover (− 28,707 km2 ± 9767 km2), we find a doubling in total areal coverage of vegetation (111% ± 13%), a quadrupling in wetlands coverage (380% ± 29%), increased meltwater (15% ± 15%), decreased bare bedrock (− 16% ± 4%) and increased coverage of fine unconsolidated sediment (4% ± 13%). We identify that land cover change is strongly associated with the difference in the number of positive degree days, especially above 6 °C between the 1980s and the present day. Contrastingly, absolute temperature increase has a negligible association with land cover change. We explain that these land cover changes represent local rapid and intense geomorphological activity that has profound consequences for land surface albedo, greenhouse gas emissions, landscape stability and sediment delivery, and biogeochemical processes.

The Arctic has been warming at double the global mean rate since the 1970s 1 .Some of the most pronounced recent warming has been across Greenland, where mean annual air temperatures between 2007 and 2012 were 3 °C warmer compared to the 1979-2000 average 2 .More extremes of temperature and precipitation are expected in the near future as Greenland's climate resilience decreases and non-linear land-climate system feedbacks develop 3 , including soil development and vegetation change, land surface albedo change, and permafrost degradation.The environmental impacts of Arctic climate change are most obviously manifest in Greenland's abundant, expanding and rapidly-evolving proglacial landscapes 4 .Specifically, Greenland's glaciers and ice caps (GICs) are shrinking, glacier-fed lakes are expanding, permafrost lakes are draining, rivers are transporting vast amounts of sediment and aggrading and widening, and vegetation cover and species diversity is expanding, largely coincident with Arctic shrubification [5][6][7][8][9][10][11][12] .
Understanding ongoing climate-landscape interactions across Greenland is crucial for modelling Arctic climate 13 , monitoring and managing water resources 14 , determining the health and livelihoods of Arctic societies 15 , and for maximising economic development prospects 16 .Specifically in Greenland, climate-landscape feedbacks include: (i) exposure of dark surfaces such as water and bedrock that have a high absorption of solar energy, i.e. a low albedo) causing increased thermal retention 17,18 ; (ii) increased deposition of organic matter, both terrestrially and within streams and eventually fjords as a consequence of warming-induced soil formation and permafrost degradation which drives further colonization of vegetation in areas where sediment redistribution is common and stream networks diverge regularly 19 ; (iii) expedited degradation of permafrost by vegetation establishment and expansion subsequently releasing substantial previously stored greenhouse gases 12,20,21 ; and (iv) increased/altered biogeochemical weathering and reactions coinciding with changes in microbial community composition driven by warming and increased moisture within the environment which affects gas (e.g.CO 2 , O 2 , N 2 O) production and sequestration rates 22,23 .The rates of change associated with these feedbacks is highly

Classifying land cover and quantifying land cover changes
We designed and implemented a rigorous pre-processing and classification workflow, as fully detailed in our Supplementary Information (SI), leveraging Google Earth Engine (GEE) computing power to analyse imagery from two time periods; Landsat-5 Thematic Mapper (TM) for the late 1980s (July to Sept, 1986 to 1989 inclusive) and Landsat-8 Operational Land Imager (OLI) for the late 2010s (July to Sept, 2016 to 2019 inclusive).Misclassification from differing illumination conditions due to topographic shading and differing solar positions during image acquisition was reduced by applying a topographic correction 34 and a pixel quality weighted mosaicking routine.The latter produced a mosaic from summer month images on a pixel-by-pixel basis using a weighted median filter based on the Near InfraRed band to reduce the inclusion of low reflectance (shadow) and over-saturated (snow, cloud edge).We implemented a semi-supervised Random Forest classification 35 to derive nine broad land cover classes; bad data (NoData/cloud/shadow), ice and snow, meltwater (e.g.rivers and sediment plumes), freshwater, coarse sediment, fine sediment, bedrock, dry tundra, and wet/dense vegetation (Fig. 1).Detailed class definitions are given in our SI.These nine classes were chosen to represent the major land cover of the biophysical environment whilst reducing redundancy (too many classes) and maintaining accuracy (distinction between classes).The overall accuracy, which we derive using a proportional area error estimation method 36 (SI), is 85% for our 1980s classification and 82% for our2010s classification.The full methodology and results of accuracy assessment, adjustment for error, and calculation of confidence intervals is outlined in the accompanying Supplementary Information.
To directly address the accuracy of changes detected we constructed a 64 class change image and correspondence matrix, for which we found change accuracy of 69%.We identified that misclassification predominantly occurs between spectrally-similar classes and error is spatially heterogeneous, so we aggregated classes with similar geophysical and geomorphological characteristics and computed spatial distribution of error via geographically weighted binomial regression (SI).
To identify spatial patterns, percentage area coverage of each of the nine broad land cover classes were spatially-aggregated on a 10 × 10 km 2 grid for each time period (SI).To understand drivers of the regional pattern of land cover changes, and due to sparse and intermittent field measurements of weather across Greenland (e.g. 37), we compared our observed large-scale patterns of land cover change to Copernicus Climate Change Service (C3S) ERA5 modelled air temperature datasets 38 .Intra-catchment processes were detected and analysed by exploiting the high spatial resolution (30 m) of our land cover classifications to identify exemplar sites of land cover change types.

Regional patterns of land cover change
Overall, we find there has been a 28,707 km 2 ± 9767 km 2 areal loss of snow/ice from the ice sheet margin and GICs, and that has predominantly been replaced by barren ground; bedrock and coarse sediment.The losses of ice volume associated with these land cover changes represent considerable contributions to global sea level rise and ocean freshening 39 and will also control river discharge and associated sediment plumes in fjords affecting coastal ecosystems (e.g. 40).As expected, we detect ice loss concentrated around the edges of present-day glacier margins, but we also here reveal that it has been especially pronounced in the north and the south-west of Greenland (Fig. 2A) (c.f. 41).Highest ice losses of − 22% (> 22 km 2 of a single 10 × 10 km 2 cell) are found around 77° N in the west, with relatively high melt rates being observed in the mid-north-west (~ − 6% areas, 70° N to 74° N) and South-east (~ − 8% areas, 62° N to 67° N) (Fig. 2A).Our geographically-weighted binomial regression analysis indicated that change from snow/ice accuracy was lower in the north west than other regions, although still averaging over 80% (SI).These patterns of ice cover loss are corroborated by Mouginot et al. 42 who reported glacier mass balance from the 1972 to 2018, finding particularly high increases in ice flux (almost 100%) and cumulative loss (~ 150 Gt) from basins in mid/north-west Greenland.Ice loss is the most marked land cover change for albedo, which is a key component of surface energy balance and hence a control on micro-climate(s).
Changes in meltwater and freshwater classes have been spatially-heterogeneous, showing a slight positive correlation with latitude on both the east and west coast (Fig. 2B), and a negative trend for freshwater ponds and lakes (Fig. 2C).Across Greenland, there has been a 15% ± 15 (10,031 km 2 ± 10,516 km 2 ) increase in meltwater area, which we interpret to represent increased to river discharge and large increases in meltwater-sediment plumes in lakes and shallow parts of fjords.Our meltwater changes (Fig. 2B) are corroborated by findings from Mernild et al. 43 , who found increasing runoff trends for 80% of Greenland's catchments between 1979 and 2014 with relatively high runoff in the south-west and north-east, as we identify in Fig. 2B.This increased meltwater river discharge impacts both local fjord/coastal waters in terms of sea surface temperature, salinity, suspended sediment and stratification, all of which influence marine ecosystems 43 .More widely, meltwater increase and consequent mobilisation of fluvial sediment can increase ocean freshening and impact Atlantic meridional overturning circulation 44 besides impacting the delicate specialist marine ecosystems of the Arctic 45,46 .
Perhaps unexpectedly, we find a decrease in freshwater area of 11% ± 22% (− 4336 km 2 ± 8421 km 2 ).However, that decreased area of freshwater is due to degrading permafrost and some episodic fluxes, such as from glacier lake drainages; for example.Finger-Higgens et al. 47 reported a 2% decrease in lake area (15% for smaller ponds) in west Greenland, and Carrivick and Quincey 48 and How et al. 49 both show changes in lakes adjacent to ice-margins.
Barren ground classes have undergone incoherent spatial changes (Fig. 2D, E, F).Notable areas with loss of bedrock and coarse/dark sediment occur in the mid/north-east and mid-west of Greenland (Fig. 2D, F) and predominantly in ice-distal locations.The net change of coarse sediment is negligible at − 0.5% ± 6% (− 1183 km 2 ± 15,455 km 2 ) but there has been a − 16% ± 4% (67,553 km 2 ± 16,682 km 2 ) loss in exposed bedrock.We attribute the latter to vegetation encroachment and growth on weathered bedrock regolith surfaces and largely via the process of shrubification due to warming temperatures 50 , as studied on Disko Island, Greenland 51 , for example, and more widely quantified herein.Indeed, losses in the sediment and bedrock classes coincide spatially with increases in both vegetation classes; dry tundra and dense/wet vegetation (Fig. 2G) and these sites are also found inside the zone of highest permafrost thaw potential 21 .Exposure time since deglaciation (and www.nature.com/scientificreports/prior to 1980) in these locations is sufficient (~ 10,000 years 52 ) for soil and rudimentary peat development to be represented in our 1980s sediment and bedrock classes, producing optimal sites for vegetation encroachment during this study period 53 .Overall, loss of barren ground represents the stabilisation of land surfaces and consequently decreased sediment fluxes, which ultimately affects mineral and nutrient export to the oceans 54,55 .Fine sediment coverage has not significantly but likely increased overall across Greenland; 4% ± 13% (4274 km 2 ± 13,239 km 2 ), but spatially corresponds with increases in meltwater (Fig. 2B and E).Fine sediment movement represents increasing sediment mobility and unstable land surfaces and is commonplace during deglaciation (cf. 55) and with thawing and drying of land.Fine sediment loss has particularly occurred in southwest Greenland in the vicinity of Kangerlussuaq ~ 67° N (Fig. 2E).In the same area, there has been an increase in exposure of bedrock and a negligible change in dry tundra vegetation, which is otherwise uncharacteristic in the south-west (Fig. 2E, F, G).These land cover changes can be explained by expansion of deflation patches, which are proceeding at rates of ~ 2.5 cm yr −156 .There is also a spatial association in south-west Greenland between this aeolian erosion of fine sediment (sand and silt) and the increases in cover of vegetation.On the east coast of Greenland, fine-grained sediment has increased at 72 o N in Jameson Land and that corresponds with reduced bedrock coverage and increased coverage of meltwater with ice-loss further inland.The exposed bedrock in Jameson Land is an extensive Late Jurassic sandstone complex 57 and increased slope instability from permafrost degradation combined with increased glacial meltwater action and deposition has likely caused the marked increase in fine-grained sediment 58 .
Both vegetation classes have in combination increased in coverage across Greenland by 111% ± 12% (87,475 km 2 ± 15,288 km 2 ).We find pronounced increases in vegetation cover across the south-west, east, and north-east (Fig. 2G).The greatest increase in coverage of dense/wet vegetation has occurred in the vicinity of Kangerlussuaq in the south-west and isolated patches in the north-east (Fig. 2H).Most of this new dense/wet vegetation coverage can be attributed to wet heath development along receding lake shores between 1995 and 2017 59 .Jørgensen et al. 29 identify similar patterns of fenland distribution in the north-east (~ 2012) as those in our contemporary classification, but studies from the mid-late twentieth century [60][61][62] indicated far less extensive dense and wetland vegetation across the north-east, as depicted in Fig. 2H.Overall, changes in both vegetation cover classes show a latitudinal gradient; a pattern of increasing change with latitude to between ~ 63 o N and 69 o N and decreased change north of this (Fig. 2G, H).Vegetation changes are especially important because they represent a stabilisation of land, introduce seasonal albedo variability and strongly affect snow retention, infiltration and overall local boundary layer (micro)climate.Moreover, vegetation development stabilises hillslopes reducing sediment and nutrient delivery to watercourses, yet simultaneously expedites permafrost degradation and introduces carbon and dissolved organic carbon and other greenhouse gases to the terrestrial 63 , fluvial 64 , oceanic 65 and atmospheric environment 66 .

Geomorphological processes
Four locations (see Fig. 1 for locations) can be chosen to exemplify landscaping processes and resultant landform and land cover changes that are representative of the spatial trends and that are permitted by the native 30 m resolution of our analysis.Figure 3A shows the southernmost source of the Qinnguata Kuussua river in south-west Greenland (Fig. 2E) where meltwater is primarily subglacially-fed and the absence of glacier surface drainage reflects low precipitation and relatively little supraglacial melt 67,68 .The proglacial river braidplain is complex in planform and composition (Fig. 3A) but is notable for the expansion in coverage of 3% or ~ 10 km 2 of the 400 km 2 selected area of fine-grained sediment, reflecting aggradation of fluvially-transported sediment.This type of braidplain expansion and altered river morphology resulting from sediment aggradation has far reaching implications.Terrestrial sediment deposition and anastomising may reflect reduced energy and so reduced nutrient conveyance marine ecosystems 69 .However, as sand reserves are depleted globally and demand increases (e.g.urban expansion, infrastructural improvements, coastal protection), Greenland may be poised to exploit this widespread phenomenon (4% increase in fine sediment across Greenland) to become a leading global aggregate exporter for sustainable economic independence 70 .
Figure 3B evidences increased sediment flux and consequent river changes around the ice-dammed lake Blåsø in north-east Greenland.There have been river channel avulsions and sedimentation across 3% or ~ 11 km 2 of the 400 km 2 area and that has most obviously occurred as delta progradation.New sediment is only found around the mouths of meltwater streams (Fig. 3B), thereby refuting attribution to a falling lake level.Increased meltwater discharge, which is ubiquitous around Greenland (Fig. 2B) and the associated fluvial mobilisation of sediment is most likely to be causing the delta progradation.Where sediment supply is abundant, i.e. regions we identify as increasing in sediment and bare ground classes largely in proglacial landscapes (Fig. 2D, E, F), the meltwater rivers transport that sediment ultimately driving widespread progradation in lakes and fjords (Fig. 3B).2][73] ).
Figure 3C shows the Saefaxi Elv river which drains Centrum Lake eastward to the ocean via Hekla Sound in King Frederick VIII Land, north-east Greenland.There has been a 4% or 16 km 2 increase in the surface area of the river over the 400 km 2 area, which can be predominantly attributed to the Greenland-wide increase in meltwater discharge 7,43 ; we find meltwater river area increases across Greenland of 15% coverage, though a large proportion of this will encompass sediment plumes into lakes and fjords, as well as some wet ice.We also note considerable increase in both vegetation classes in this area (Fig. 3C), and evidences the spatial association of vegetation encroachment and growth with increased availability of freshwater sourced from melting GICs and degrading permafrost 74 .
We find a near quadrupling (increase of 380% ± 29%, 30,295 km 2 ± 6805 km 2 ) in the coverage of wet/dense vegetation across Greenland.Wetland expansion on Shannon Island, east Greenland (Fig. 3D) is an exemplar of widespread wetland development and expansion in east and north-east Greenland (Fig. 2H).Jørgensen et al. 29 identified wetlands and fen in this region and proposed that the wetlands would continue to expand with permafrost degradation.These wetlands are likely minerotrophic, principally deriving nutrients from glacial meltwater, and constitute considerable sources of CH 4 in polar regions 75 .Locations of wetland formation and growth in north-east Greenland identified here present key sites for research to explore the poorly understood www.nature.com/scientificreports/process of wetland, and peat, initiation and lateral expansion, a field at the forefront of ongoing arctic biospheric research 76,77 .Arctic wetland development has a complex and poorly understood biogeochemical connotation for greenhouse gas emission versus drawdown; Arctic wetlands are effective methane emitters 78 , yet dry and moist arctic tundra soils in north east Greenland have been measured to draw down methane 29 .Expansion of vegetation and especially in wetland areas indicates but also exacerbates permafrost thaw, active layer thickening and thus emissions of greenhouse gases previously stored in these Arctic soils 79 .

Association of land cover change with climate change
We conducted a number of spatial regression analyses to determine how climate change is associated with the land cover changes.The difference in the average (per decade) number of days per year above various significant positive degree temperature thresholds was more informative than absolute magnitude of temperature change because of the confounding influences of temperature and precipitation.Consider a three degree warming that remains below zero versus a three degree warming that crosses from below to above zero thereby leading to ice melt and permafrost thaw, for example.Comparing our land cover change grids to our Difference in Degree Days Above Temperature grids (between 1980 and 2010s) on a grid cell by grid cell basis, we found lowest R 2 and Akaike information criterion values (SI) for the two water and the fine sediment classes.We find the highest R 2 and Akaike information criterion correlations for the vegetation, bedrock and ice classes, suggesting a close association between these types of land cover change and climate change across Greenland.Specifically and most interestingly, the highest R 2 and Akaike information criterion values for all classes are found for the difference in the number of days above 6 °C.This leads us to suggest that future land cover change is therefore likely to be marked and accelerated in regions where the number of days per year above this 6 °C threshold are increasing most rapidly; crucially that is not simply where temperatures are set to rise most rapidly.Table SI 12 contains all geographically weighted regression results.

A model of changing land cover across Greenland
We summarise both regional inter-class changes and intra-catchment geomorphological processes that we have measured (Figs. 2, 4A) within a conceptual model (Fig. 4A, B) that is based upon 64 class change matrices.Each of those 64 class changes is a calculation of every 'from-to' change type between classes (Fig. 4C), and considering both relative change within each matrix as well as relative change per class between matrices (SI).We highlight that there is a notable difference in the latitudinal pattern of the predominant land cover inter-class changes per latitudinal band between the east and west of Greenland (Figs. 2, 4B).The size of proglacial areas is far larger on the west coast than in the east, and the topographic relief is generally more subdued in the west compared to the alpine east and so accommodating of regolith and tundra vegetation establishment.Moreover, the two coasts have experienced different climatic change 37 , in terms of temperature and precipitation; the east is influenced by variability in North Atlantic Oscillation and the Atlantic Meridional Overturning Circulation.
Our model, and the relationship between land cover changes and the difference in degree days above 6 °C, illustrates that higher air temperatures are producing more suitable conditions for vegetation expansion and growth 50 .Vegetation cover is perhaps the most profound land cover change we have quantified because it has a considerable feedback effect on the climate system 80 .For example, the newly-established and enlarging wetlands and fens in phase 6 regions of the north-east of Greenland are associated with considerable methanogenesis 81 , but that is potentially partially offset by long-term carbon storage via peat accretion 82 .Field observations combined with mapping indicate widespread shrubification in Greenland, whereby the "greening" of the Arctic is largely manifest in the proliferation of shrub tundra species 83 .
We identify substantial vegetation areal expansion and growth in regions of critical permafrost thaw (phases 5 and 6, Fig. 4).That expansion can be interpreted as a manifestation of increased air temperature warming ground surfaces and driving permafrost thaw 84 .Permafrost thaw impacts a suite of human and economic activity; Greenland's infrastructure has been forced onto challenging ground underlain by thick poorly consolidated glacigenic sediment and permafrost 21 with economic development, localised population growth and increased tourism that all have environmental consequences (c.f. 85).Our analysis can therefore be interpreted to reveal a positive feedback between permafrost degradation and vegetation encroachment that will likely increase exponentially, threating both local communities and infrastructure 15 .These land cover feedbacks, especially for landscape parts in phases 5 through 7 (Fig. 4), therefore warrant policymakers' attention.
Our study has identified, quantified and summarised regional, inter-catchment and intra-catchment land cover changes across Greenland between the late 1980s and the late 2010s.This spatiotemporal coverage has enabled proposal of a conceptual model of land cover phase change and the association of that land cover change with climate.Specifically, we find that where the difference in degree days is > 6 °C, then we find increased icesheet and glacier mass loss, increased meltwater production and increased sediment mobility.Overall, the rapid changes in land cover across Greenland requires repeated monitoring to identify change trajectories and to refine understanding, which will directly benefit human and economic development.For the geosciences, feedbacks of land cover with climate must be considered when modelling and projecting climate, the Greenland Ice Sheet, GIC evolution, meltwater runoff, sediment yield and habitat and ecological community changes.

Methods
Here we briefly outline the main methods for producing our image mosaics and classification., it being freely available, and its relatively high 30 m resolution [87][88][89] .TOA products were chosen over atmospherically corrected surface reflectance (SR) as the US Geological Survey (USGS) 90 state the efficacy of surface reflectance correction is reduced in: (i) highly snow-covered regions, (ii) regions with low-sun angles, (iii) coastal regions where land area is small relative to adjacent ocean, (iv) areas with high cloud cover conditions, and (v) high latitude areas over 65° North 90 .As Greenland's peripheral ice-free areas meet most/all of these criteria, SR products should not be used for analysis.Summer month images (July-September) are selected for the years 1986 to 1989 inclusive for the late eighties and 2016-2019 inclusive for the contemporary mosaic.We produced multi-month multi-year mosaics as: (i) the temporal resolution is limited by the 16-day return period for Landsat imaging, (ii) cloud cover and orographic shadowing is particularly problematic in certain regions of Greenland 91 , and (iii) there is poor Landsat coverage for the south, south west, and north west of Greenland during the late eighties 92 .Full workflow is outlined in (Figure SI.1).

Image classification
In the first instance, the classification was conducted using an unsupervised, k-means clustering algorithm executed in GEE.This method removes the users potential to target predefined purpose-driven classes which may be difficult to identify or that are presumed to be in abundance when relatively spatially sparse 93,94 .The unsupervised approach is particularly useful for large, national/regional scale spatial analysis where it is unreasonable for the user to cover the entire area for validation and therefore only classes which can be well distinguished and are spectrally distinct are produced.The initial classification workflow is presented in Figure SI.4.Once a classification is produced with K number of classes, these are then aggregated into 9 land cover classes.These final 9 land cover classes were chosen as they are deemed distinct and identifiable from the imagery, present across the entire periphery of Greenland, relevant to this research and enough as to reduce the production of redundant classes.These classes and their descriptions are shown in Table SI. 1.

Classification accuracy assessment
Accuracy of the contemporary classification was assessed using four band (RGB,NIR) 3 m resolution contemporary PlanetScope basemap mosaic imagery from the third quarter of 2019 (July-September) as reference 95 and a Sentinel-SA image mosaic for 2019 summer season produced in Google Earth Engine.Due to a lack of supplementary satellite imagery for the 1980s, our Landsat mosaic was used for reference alongside 2 m resolution greyscale georeferenced aerial orthophotographs, made available by the AeroDEM project 96 .We report an overall accuracy for our classifications to one significant figure only; 85% for our 1980s classification and 82% for the 2010s classification.Below we outline in detail our accuracy assessment procedure using the contemporary classification.This same method was also applied to our eighties classification.Accuracy assessment was conducted for six 100 km × 100 km areas south of 76° N shown in Figure SI.2.Plan-etScope imagery is not available north of 76° N. Stratified random sampling based on total occurrence of classes for the entirety of Greenland was used to generate 500 accuracy assessment points (AAP) distributed between the 8 land cover classes at each location.Results were then amalgamated resulting in 3000 AAP over total area of 60,000 km 2 .We removed the ice sheet from analysis as this skewed the stratified AAP production and is largely static over the study period away from the margins.Snow and ice accuracy was therefore assessed using GICs.A correspondence (AKA: error, confusion) matrix was produced using the 3000 AAP and is presented both absolutely and as an expression of proportional area considering the associated error within each class, following Olofsson et al. 36 .The overall kappa coefficient is 0.781, and overall accuracy (OA) was presented as 82%.The basic kappa coefficient however adjusts for 'random allocation agreement' and the validity of the adjustment has been questioned 98,99 .We therefore also used a stratified estimation of error following Olofsson et al. 36 to mitigate against potential inherent measurement bias and to estimate standardised errors with quantified uncertainty based on sampling variability.The cell-count derived correspondence matrix is presented in Table SI.2.

Field validation
In order to further validate our contemporary classification we collected field validation at a spatially distributed array of 89 sites as mapped and detailed in Figure SI.6.Photographs of each land cover type are presented in Figure SI.7.We found our field validation points to have an overall kappa coefficient of 0.81, in line with the accuracy we reported from remote validation.

Classification change accuracy assessment
To confirm that the accuracy of the individual classifications is reflected in the changes we measure, we directly quantify the accuracy of changes.Using Google Earth Engine (GEE) a 64 class change map for the entirety of Greenland was produced (available as supplementary dataset).See Table SI.5 for legend describing the 64 change class designations.
A further 3000 AAP were randomly generated, stratified by the occurrence of each of the 64 classes in Google Earth Engine.The minimum number of points for each class was set at 5, as some classes has especially low areas which equated to less than 1 point when stratified by relative areas (e.g.class 56: Dense vegetation to Snow/ ice).The accuracy of these changes was then assessed using the same reference data as was used for individual class accuracy assessment, i.e.Sentinel-2A mosaic and Planet imagery for contemporary, Landsat mosaic and 2 m AeroDEM orthophotographs for the 1980s.The overall Kappa statistic of accuracy for the change image correspondence matrix 0.69.Whilst conducting change class accuracy assessment it was evident that accuracy was not spatially homogenous, and misclassifications were predominantly between classes with largely similar spectral signatures and associated geomorphological process attributions.We therefore explored the spatial distribution of error following class aggregation as outlined in Table SI.6.The specifics and findings of spatially distributed error estimation are outlined in our SI.

Associations with climate
We conducted a number of regression analyses to understand the degree to which Greenland's climate change has influenced the land cover changes we observe over the thirty year study period.Decadal and annual average summer temperature change from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-5 monthly average dataset for summer months was correlated with the change classes, but Ordinary Least Squared (OLS) correlations were all below 0.1 and indicated a non-parametric relationship.We speculated that average increase in temperature is not a useful an indicator of land cover in Greenland because changes above 0 °C will produce ice melt and permafrost thawing, whereas a change (of equal magnitude) below zero will not.We therefore calculated the difference in the average number of days in a year above three significant positive degrees temperature thresholds (0 °C, 3 °C and 6 °C) for the late 1980s and late 2010s which we term Difference

Figure 1 .
Figure 1.Landsat-8-derived land cover classification for the 2010s.The coloured circles and site names indicate where we highlight analysis of catchment-scale processes in Fig. 3, and the black squares with letters denote five locations within which are our 89 field validation points detailed in Fig. SI.6.

Figure 2 .
Figure 2. Percentage of total area change within 10 km × 10 km grid cells for each land cover class and plots of variation in average land-cover change with latitude for both east and west Greenland and 95% Confidence intervals.Inset K Kangerlussuaq, JL Jameson Land.

Figure 3 .
Figure 3. Examples of localised land cover change demonstrating geomorphological processes driving environmental evolution as seen using changes in 30 m resolution land cover data.(A) Proglacial sedimentation (Qinnguata Kuussua river), (B) delta progradation (Blåsø lake), (C) increased meltwater stream size and discharge (Saefaxi Elv), and (D).vegetation and wetland growth (Shannon Island).SS = Suspended Sediment.Locations for sites in panels are shown in Fig. 1.
https://doi.org/10.1038/s41598-024-52124-1www.nature.com/scientificreports/Image preparation Landsat-8 Operational Land Imager (OLI) Top of Atmosphere (TOA) imagery is used for the contemporary classification, and Landsat-5 Thematic Mapper TM TOA imagery is used for the late 1980s classification.Landsat imagery is the most widely used data type for land cover mapping and classification largely due to it having the longest satellite image record spanning nearly 50 years

Figure 4 .
Figure 4. Land cover change model 1980s to 2010s for Greenland, illustrating numbered predominant landcover inter-class changes on a common trajectory reflecting increasing landscape stability with time; i.e. phases (A).Panel (B) shows degree day difference air temperature (DDDAT) mapped and with east and west coast latitudinal transects of numbered predominant landcover inter-class change in each latitudinal band.Numbered inter-class changes are detailed in the table (C). https://doi.org/10.1038/s41598-024-52124-1