Glacier Retreat in Iceland Mapped from Space: Time Series Analysis of Geodata from 1941 to 2018

In recent decades, glaciers outside Greenland and Antarctica have shown increasingly rapid rates of mass loss and retreat of the ice front, which is associated with climatic and oceanic warming. Due to their maritime location, Icelandic glaciers are sensitive to short-term climate fluctuations and have shown rapid rates of retreat and mass loss over the last decade. In this study, historical maps (1941–1949) of the US Army Map Service (AMS series C762) and optical satellite imagery (Landsat 1, Landsat 5, Landsat 7, Landsat 8, and Sentinel-2) are used to study the Langjökull, Hofsjökull and Vatnajökull ice caps. By the help of the Normalized Difference Snow Index (NDSI), the glacier terminus fluctuations of the ice caps from 1973 to 2018 and the Equilibrium Line Altitude (ELA) from 1973 to 2018 are analyzed. The results are compared with climate data, especially with mean summer temperatures and winter precipitation. Due to the negative temperature gradient with increasing altitude, bivariate histograms are generated, showing the glaciated area per altitude zone and time, and providing a prediction of the future development until 2050 and beyond. The results indicate that Langjökull, Hofsjökull and Vatnajökull are retreating and advancing over the study period in correlation with the mean summer temperature, with a steady decrease over time being the clearest and most significant trend. The lower parts of the glaciers, thus, will probably disappear during the next decades. This behaviour is also evident by an exceptional increase of the ELA observed on all three glaciers, which leads to a reduction of the accumulation zone.


Introduction
Glaciers and ice sheets cover about 11% of the Earth's surface. They serve as drinking water reservoirs and are the basis for hydroelectric power, especially in Iceland. According to the Intergovernmental Panel on Climate Change (IPCC) (Intergovernmental Panel 2021), glaciers are losing mass in response to the clear warming of the atmosphere and the oceans. This loss of mass has already contributed to the global rise of the mean sea level and will probably lead to a sea level rise of 1 m by the end of this century (Church et al. 2013). The global melting of glaciers can cause hazards such as the release of ice-dammed water (the so-called Jökulhlaups as observed in New Zealand and Greenland in 2007 as well as in Iceland in the year 2010) and ice avalanches (Carrivick 2011). Recent studies on the mass balance of glaciers and their ice front positions outside Greenland and Antarctica show that the rates of mass loss and retreat of the ice front have accelerated since the 1970s (Carr et al. 2014;Cogley 2009;Gardner et al. 2013;Jacob et al. 2012;Marzeion et al. 2012). Due to their maritime location, Iceland's glaciers (cf. Fig. 1), in particular, are highly sensitive climate indicators and react immediately to shortterm (annual to decadal) climate fluctuations (Bradwell et al. 2013;Sigursson et al. 2007). Studies dating from the last ten years have shown that Icelandic glaciers have experienced rapid rates of retreat of the ice front and mass loss Bradwell et al. 2013;Sigursson et al. 2007;Hannesdóttir et al. 2015). The mass balance of glaciers is a key variable that needs to be monitored for a better understanding of the interrelation between glaciers and climate. It is traditionally monitored using glaciological methods based on point, pit, and probe measurements. With the increasing availability of remote sensing data, the glacier mass balance can also be assessed by identifying the snowline on images taken at the end of the ablation season for each year, which can be considered as the equilibrium line (EL) for temperate zone glaciers (Paterson 1994).
At the end of a mass balance year, just before the first snowfall, the boundary between accumulation and ablation zones becomes visible. The EL forms a balanced line between accumulation and ablation. Therefore, the altitude of the EL (ELA) is an important indicator for the local climate as well as an important feature of the glacier mass balance. Depending on the weather conditions, either the firn line (minimum firn height, marking the altitude of the long-term EL) or the transient snow line (minimum snow height, marking the ELA at the end of an ablation season) will be identified (Jaenicke et al. 2006). The end-of-season snowline height (SLA) is an indicator of the ELA of a glacier and is often used to derive its mass balance. When the mass balance of a glacier is negative, i.e., ablation is higher than accumulation, the ELA increases. If the mass balance is positive, the ELA moves to lower altitudes (Kuhn 1981). In general, the growth and decay of glaciers assessed by the surface mass balance depend on climatic conditions. The analysis of glacier parameters such as area, surface zones and mass balance of Icelandic glaciers is, therefore, crucial to estimate the current period of atmospheric warming and the associated glacier retreat.
An early study of Iceland's glaciers initially focused on glacier variations from 1930 to 1995 derived from the database of the Icelandic Glaciological Society published in 1998, which contains about 60 glaciers and glacier tongues (Sigursson 1998). Sigursson et al. (2007) gives a brief overview of available data on terminus fluctuations of Icelandic glaciers, since regular monitoring started in the 1930s, and shows examples of how they reflect temporal variations in the climate of the island. During the twentieth century, decadal temperature variations in Iceland had a significant influence on the variations of the glacier termini. A correlation between variations in the mass balance of the glacier and changes in summer temperatures was observed (Sigursson et al. 2007). For more than one century, the World Glacier Monitoring Service (WGMS) World Glacier Inventory (1988) and its predecessor organizations have been collecting and disseminating standardized data on glacier fluctuations. In close cooperation with the U.S. National Snow and Ice Data Center (NSIDC) (2021) and the WGMS, the initiative Global Land Ice Measurements from Space (GLIMS) (Raup et al. 2007) was established to continue this inventory task with space-borne sensors. They have mapped all glaciers in Iceland in the period 1999-2004 using satellite, map and aerial image data. The outlines of all glaciers were recorded, digitized on screen and added to a global glacier database. While in the past, optical remote sensing data were used exclusively to delineate the equilibrium line using spectral indices or band ratios, a synergy of optical and radar remote sensing data was recently used to monitor and evaluate the mass balance of the Mýrdalsjökull ice cap in Iceland (Jaenicke et al. 2006  . Using ice-front measurements, photographic and geomorphological evidence, the recent retreat of Virkisjökull-Falljökull, south-eastern Iceland, was examined (Bradwell et al. 2013). Hannesdóttir used multiple digital elevation models from glacial geomorphological features (including lateral and terminal moraines, trim lines and glacier erratics), historical photographs, maps, aerial images, DGPS measurements, and a LiDAR survey to analyze area, volume and mass changes of the southeastern Vatnajökull ice cap from the Little Ice Age maximum in the late nineteenth century till 2010 (Hannesdóttir et al. 2015). The Icelandic Meteorological Service (Icelandic Met Office 2021) also maps the changes of the Hofsjökull and Drangajökull glaciers twice a year, in spring and autumn.
Various methods have been developed to identify glacier dimensions from imagery. The detection of a clean glacier surface is based on the spectral uniqueness of glacier ice in the visible and infrared part of the electromagnetic (EM) spectrum. According to Rees (2006), snow and ice are characterized by: (1) high reflectivity (albedo) in the visible wavelengths (0.4-0.7 μm); (2) medium reflectivity in the near-infrared (0.8-2.5 μm); (3) low reflectivity and high emissivity in the thermal infrared (2.5-14 μm); and (4) low absorption and high dispersion in the microwave signals. Single band ratios provided by optical sensor systems Kääb et al. 2002;Bayr et al. 1994) such as the Normalized Difference Snow Index (NDSI) are frequently applied means to detect and classify ice cover (Hall et al. 1995). The NDSI has proven to be a robust method (Hall et al. 1995;Racoviteanu et al. 2008;Sidjak 1999). In a recent study (Nagler et al. 2016), the advancements for snowmelt monitoring by means of Sentinel-1 synthetic aperture radar (SAR) were reported. While a major advantage of the C-band SAR for snow mapping is its ability to penetrate clouds, its sensitivity to dry snow cover is rather limited. Therefore, a combination of optical sensors and SAR is recommended (Nagler et al. 2016). Nevertheless, Heilig et al. (2019) recently published an approach to derive the wet snow zone from time series of space-borne SAR sensors exclusively. They identify the wet snow by its specular scattering, which reduces the measured backscattering intensity to a minimum. The outlines of the wet snow zone at the end of the ablation period consequently define the EL. The intersection with a DEM leads to the ELA. Therefore, SAR data could be useful in the future. However, today's archives are often still too sparse and cover a too short period of time to derive climate-relevant information.
Long and continuous records of glacier fluctuations are an integral part of monitoring glaciers at a regional and global level. Unfortunately, sufficiently detailed records may be missing or at least incomplete in many places all over the world. Previous studies of Iceland's glaciers have mainly focused on monitoring special glacier parameters and/or individual glaciers (Bradwell et al. 2013;Hannesdóttir et al. 2015;Björnsson et al. 2013;Jóhannesson et al. 2013). Only a limited amount of extended time series analyses have been carried out with their focus on glacier termini fluctuations and climate parameters. Hence, a prolonged time-series analysis of multiple glaciers in Iceland with the above-named focus is beneficial to assess the response of a glacier to climate change.
Therefore, the aim of the present work is to analyze the glacier fluctuations of Langjökull, Hofsjökull and Vatnajökull over a long period from 1941 to 2018 and their possible relationship to climate parameters, as well as to determine the movement of the ELA between 1973 and 2018 since remotely sensed image data have been available. We are well aware, that the height information needed for a distinct mass balance is drawn from a potentially outdated DEM. But, the frequent acquisition of high-resolution DEMs is very expensive, whereas the space-borne imaging systems used in this study provide a weekly to monthly coverage with high spatial resolution at no cost. Hence, one question will be, whether the study of the freely available satellite images-and thus, 2D information-might lead to comparable results as the studies using cost-intensive height models and in situ data gathered during field campaigns.

Iceland's Glaciers
Iceland is a partly glaciated, volcanically active island in the North Atlantic, and lies just over 300 km east of Greenland, 1000 km west of Norway and 800 km northwest of Scotland. Iceland has an area of 103,000 km 2 , with glaciers covering about 11% of its land area. Vatnajökull, Langjökull and Hofsjökull represent the three largest ice caps of Iceland in descending order (cf. Fig. 1). The Vatnajökull ice cap partially overlaps an active volcanic zone on its western half, whereas Langjökull and Hofsjökull lie completely within an active volcanic zone (Björnsson et al. 2003). The Grímsvötn volcanic system on western Vatnajökull is one of the volcanoes in Iceland with the highest eruption frequency and the highest volcanic productivity (Thordarson and Larsen 2007).
Glaciers generally can be classified into surge-type and non-surge-type glaciers with regard to the reaction of the terminus to variations in mass balance. The terminus of a non-surge-type glacier will react predominantly to changes in the mass balance (Sigursson et al. 2007). Non-surging areas are in the southeastern part and the northwestern flank of Vatnajökull, the northern, western and eastern outlets of Langjökull as well as the northern part of Hofsjökull ) (cf. Fig. 1). On the contrary, the terminus of a surge-type glacier retreats most of the time regardless of climate and mass balance. It only advances occasionally at a more or less regular interval of time (Björnsson et al. 2003). The reaction of the terminus of a surge-type composite glacier, coming from source areas with variable characteristics, can be influenced by factors other than the climate. Corresponding factors may be geothermal or volcanic activity, debris on the glacier surface and calving in a lake or the ocean. Surge-type glaciers occur almost exclusively as outlets of the major ice caps. However, several surge-type glaciers also exist outside of the major ice caps: three of them can be found in the mountains of north-central Iceland, for instance. On the Vatnajökull glacier, approximately 75% of its outlets are surge-affected. Surges tend to occur in the southern outlets of Langjökull and especially in the flat ablation areas of eastern Hofsjökull (cf. Fig. 1).
However, previous studies report on how the mass balance of glaciers can also control the frequency of surges (Harrison and Post 2003). A change in a glacier's climate and mass balance could prevent glaciers that have surged in the past from accumulating mass to initiate further surges and thus, bring them out of the surge cycle (Hodgkins et al. 1999;Hansen 2003;Malecki et al. 2012). Striberger et al. (2011) predict that the surging of Eyjabakkajökull (eastern Iceland) and similar surge-type outlet glaciers will no longer occur frequently in comparable areas due to reduced net mass inflow rates and the extension of the quiescent phase (the period of slow ow and/or stagnation for tens to hundreds of years between surges) in a warming climate, and that Eyjabakkajökull, in particular, may eventually even fail to return to its active phase and possibly switch to a nonsurging mode. Sevestre and Benn (2015) suggest that the existence of climatic envelopes that favour surging implies that glaciers can change from "normal" to surge glaciers and vice versa under cooling or warming climatic conditions. According to Björnsson and Pálsson (2008), Iceland's glaciers are classified as "warm-based" or "temperate" and react actively to climate fluctuations. As temperate glaciers exhibit surging behaviours, it is assumed that the highest density of surge glaciers occurs in a relatively narrow climatic band limited by mean annual temperatures of about − 10° to 0° and mean annual precipitation of about 200-2000 mm (Sevestre and Benn 2015). With its maritime cold temperate to low arctic climate and numerous temperate glaciers, Iceland lies largely within this climate band. Surge-type glaciers in Iceland are both hard and soft bedded, overlay a variety of volcanic substrates and are located in regions that are exposed to a variety of climatic conditions (Björnsson et al. 2003).
Climate-related glacier changes in Iceland are driven by the elements of climate, e.g., temperature and precipitation (Meier 1965;Bradwell et al. 2013), although these depend on the topographic altitude and terrain (Einarsson 1984). In Iceland, several non-climatic factors may in some cases also influence the mass balance. Such factors can be sub-glacial geothermal or volcanic activity, debris on the glacier surface, volcanic air fall (tephra) and iceberg calving into water bodies (Sigursson et al. 2007). While both summer temperatures and winter precipitation vary, the long-term trend indicates an increase in both climate parameters (cf. Fig. 2). In recent years, from 1990 onwards, except for Vatnajökull in the last decade, an overall negative trend of glaciated area is detectable (cf. Figs. 3, 4).

Historical Maps
Langjökull, Hofsjökull, and Vatnajökull are recorded on a total of 45 historical maps of the US Army Map Service (AMS series C762). The Danish General Staff (later Danish Geodetic Survey) has begun plane table mapping of Iceland including its glaciers in 1902 to establish a new triangulation of the Reykjanes peninsula in the southwest (Bövarsson 1996). Based on this early field survey, distance measurements, and aerial photographs by the US Army Map Service (acquired between 1945 and 1946), the U.S. Department of Defense published several map series at different scales of Iceland, including the C762 series used in this study with full coverage on a scale of 1:50,000 in cooperation with the Iceland Geodetic Survey (Landmaelingar Íslands). The maps are given in Transverse Mercator projection. The horizontal date refers to the Reykjavik astronomical station at 64 • 8 ′ 31.88 ′′ northern of the equator and 21 • 55 ′ 51.15 ′′ western of Greenwich. There are two coordinate grids with a spacing of 1 km plotted on the maps: the Icelandic grid (as described above) drawn in brown and the Universal Transverse Mercator (UTM) grid (zone 27 N or 28 N) drawn in black. The maps were georeferenced via the coordinate grid and then digitized in preparation of this study.

Satellite Images
A total of 58 multispectral remote sensing images, including 48 Landsat (Landsat 2021) and 10 Sentinel-2 (Sentinel 2018) images, dating between 1973 and 2018, were obtained. Landsat imagery was acquired by Landsat 1 Multispectral Scanner (MSS), Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper (ETM+) and Landsat 8 Optical Land Imager (OLI). Only images acquired towards the end of the ablation season (June-September) were used due to the minimum snow cover. Additionally, the area of interest had to be cloud-free. The Visible Green and Short Wave Infrared (SWIR) bands were used to calculate the Normalized Difference Snow Index (NDSI). In the absence of SWIR channels, the near-infrared (NIR) band was used for Landsat 1 (cf. Table 2).
Sentinel-2 data was available in Level-2A and is, thus, delivered as bottom of atmosphere (BOA) reflectance imagery. Landsat data are provided as Level-1 data product. The digital numbers (DN) of the Landsat images are therefore converted to top of atmosphere reflectance (TOA) first. Furthermore, the image-based atmospheric correction method DOS1 (Dark Object Subtraction) was applied. The remaining pre-processing included the projection to the same earth-fixed coordinate system UTM, which is essential for the inter-annual comparison. The pixel size of the Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI imagery corresponds to 30 m in this study. By contrast, Landsat 1 MSS is provided in a spatial grid of 60 m and the Sentinel-2 scenes possess a spatial sampling of 10 m in the Green and 20 m in the short-wave infrared (SWIR) channel. To equalize the spatial resolution of all sensors, the Landsat MSS scenes are resampled to a uniform cell size of 30 m using the nearest neighbour resampling method (United States Geological 1980). The Sentinel-2 scenes are subsampled to 30 m pixel size after low-pass filtering. Landsat 7 ETM+ data are also included in this study dating from before and after May 31, 2003, i.e., the Scan Line Correction (SLC) outage. The SLC failure produces a zig-zag pattern of void pixels perpendicular to the satellite ground track, thus an amount of around 22% of the scene is lost. To bridge the SLC-off data gaps, the data can either be merged from multiple exposures or interpolated from neighbouring pixels (United States Geological 2003). Since Iceland generally is covered by only a few cloudless scenes per year, any available image information, though perforated, is valuable. The satellite imagery for our study hence were mosaicked from multiple acquisitions within one month.
Due to the clear distinction of snow and ice from rock and soil by their spectral reflectance properties, glaciated surfaces can be extracted quite reliably using optical satellite imagery except for debris or tephra covered zones (Paul et al. 2004).
For instance, open ice shows a high degree of reflection in the visible range (0.4-0.7 μm) and medium to low reflectance in the near-and short-wave infrared (0.7-2.5 μm) (Renaud et al. 2009). In this work, the Normalized Difference Snow Index (NDSI) is used to derive the glaciated area, which is the difference between the Green band and SWIR band in relation to the sum of both (cf. Eq. 1). Since different satellites and sensors are used, the input channels of the NDSI will be adjusted accordingly (cf. Table 2). For Landsat 1, thus, the SWIR must be replaced by NIR for the lack of short wave infrared bands. The alternative NDSI (aNDSI) composes of the difference of Green and NIR (nearinfrared) divided by their sum (cf. Eq. 1). The values of the resulting output images range in-between [−1, +1] , i.e., a closed value range is produced for both snow indices by According to the typical spectra of snow and rock derived from Appendix Fig. 8 and listed in Table 1, the NDSI allows for a clear distinction of both land cover classes: 0.81 for snow vs. − 0.11 for rock. The class distinction based on the aNDSI (Eq. 1) only shows a difference of 0.24 compared to 0.92 in the case of the NDSI. With respect to the end of the summer, the border area between glacier and rock is dominated by firn or upper glacier ice instead of snow. For these two classes, the distinction capability rises to 0.32 in the case of rock and firn and even 0.58 for glacier ice and firn. Only dirty glacier ice, e.g., a tephra covered glacier, has lower values similar to snow (cf. Appendix Fig. 8). In summary, the distinction of glaciated areas from bare rock is possible with aNDSI, though the known NDSI would promise clearer results. A threshold-based method was used to distinguish snow-or ice-covered from snow-or ice-free areas (cf. of aNDSI and the traditional NDSI as well as the channel's bandwidths slightly differ from sensor to sensor (cf .  Table 1), the threshold has to be adapted as listed in Table 2.

Digital Elevation Models
The ArcticDEM provides the terrain height for the complete landmass of Iceland. ArcticDEM is an initiative of the National Science Foundation (NSF) and the National Geospatial-Intelligence Agency (NGA) for the automatic generation of high-resolution, digital surface models (DSM) of the Arctic. In this study, it is utilized to analyze the glaciated area by elevation zone as well as the ELA. The data is generated based on optical DigitalGlobe track and cross-track images with high resolution (approximately 0.5 m ground sampling distance (GSD)) (Porter et al. 2018). The elevation data obtained corresponds to a GSD of 2 m and was created in 2018. The acquisition technique (e.g., LiDAR, SAR, terrestrial photogrammetry) is known to have a major influence on the estimated elevation. The height of LiDAR, for instance, reproduces the height of the scattering center in the visible or near-infrared bands (depending on the utilized sensor) for single-point measurements. SAR interferometry models the height of the phase center in the respective microwave band in relation to the neighbouring pixels. Photogrammetry delivers the height of visible objects that can be detected and assigned to each other in separated images taken in the visible range of the EM spectrum. The DEMs derived from these methods can differ considerably (Paul et al. 2017), which complicates any multi-sensor analysis. As there is no consistent DEM timeline available for Iceland, we focus on DEM-independent detection methods based on georeferenced satellite imagery. Only the retrieved glacier outlines and ELs are finally interpreted by the help of the ArcticDEM described above being well aware that glaciers typically show certain temporal variations in height. As the highest resolution in vertical dimension required for our study is about 100 m (cf. Sect. 4.2), the influence of changing heights is more or less negligible, though a minimal impact cannot be denied.

Meteorological Data
The Icelandic Meteorological Office provides climatic data from various measuring stations and periods. The Stykkishólmur measuring station has the longest and most uniform measurement records in Iceland, dating back to 1845. It is a manned synoptic station, located in the southwest of Iceland ( 65 • 4 ′ 27 ′′ N, 22 • 44 ′ 2 ′′ W ) at 13 m above sea-level [27]. The climate tables provide information on monthly and annual rainfall, minimum and maximum temperatures, wind speed, cloud cover, barometric pressure, and several statistical averages. The data obtained were used to calculate the mean summer air temperatures (June-September), the total amount of precipitation in winter (October-April) as well as other climatic parameters (cf. Fig. 6).

Methodology
The investigation period is defined from 1941 until 2018. In principle, solar radiation and air temperature in summer, as well as, precipitation in winter are known to be important indicators for the interpretation of glacier changes in Iceland, whereby the mean summer temperature is considered the most essential parameter (Sigursson et al. 2007). Hence, the intra-annual study period is set to June-September. The analysis starts with the detection of glaciated areas in order to generate multi-temporal glacier masks. These are vectorized to derive the glacier outline for the years under study. From the cleaned outlines the glaciated area can be calculated for the different years and then plotted together. To quantitatively assess the relationship between glacier termini fluctuations and climate parameters, the glaciated area of each glacier per year is calculated. Since the elevation is assumed to play a key role as well, multivariate histograms are produced, which reflect the glaciated area by time and elevation. The quantitative assessment of the height dependency is performed by the correlation with several climate parameters. A regression is applied in order to describe the glacier change analytically. The prediction by three different models tries to forecast the glacier retreat in the coming decades. Finally, the rise of the ELA is considered to confirm the conclusion drawn from the analyses described above.

Glaciated Area
To further elaborate the glacier outlines, the glacier masks derived from the aNDSI images are converted to vector data. This is followed by the manual post-processing of the Skeið arárjökull glacier tongue (southern Vatnajökull), which is hard to detect automatically due to Tephra deposits. Accordingly, all files are checked visually before further processing and edited manually by on-screen-delineation (Hall et al. 1992;Williams et al. 1997;Paul et al. 2002), if necessary. The vectorized outlines are superimposed on the satellite image optical scenes of each year (cf. Fig. 3). The contrast of the respective satellite images is increased by a local histogram stretch to distinguish the glacier boundary from the bare soil. Subsequently, the surface areas of the individual segments are calculated. Non-glacier segments are eliminated by removing areas smaller than 31 km 2 . This threshold was chosen because Hvítárvatn, located in the Icelandic highlands (east of Langjökull) and being the largest lake in the study area, covers an area of 30 km 2 . The variation of the total glaciated area is plotted in Fig. 4 for each glacier under study. Decadal linear trends are estimated as well as a quadratic regression over the whole time frame.

Glaciated Area by Elevation
It is reported from studies on glaciers in the French Alps that the altitude and aspect of glaciers have an influence on their behaviour in the face of climate change (Rabatel et al. 2013). However, the mathematical relationships are not that easy to assess. Due to the availability of the ArcticDEM, we decided to additionally take the height into account by creating a bivariate histogram of glaciated area per elevation zone and time. The steps in ellipsoidal height are set to 100 m for stability reasons. Hence, temporal height changes play a minor to a negligible role in this study. The temporal sampling is adapted to the availability of satellite images as described in Sect. 4.1. This kind of presentation illustrates the dependence of the glacier retreat from the elevation zone over time (cf. Fig. 5).

Correlation with Climate Parameters
The multi-temporal extents, i.e., the glaciated surface area per year and elevation zone of the three glaciers under study (derived according to Sect. 4.1) are then correlated with time, the respective climate parameters like mean temperature, total precipitation, air pressure, cloud cover, and wind speed for the ablation (June-September) and the accumulation (October-April) season separately as well as with the annual total sunshine hours (cf. Sect. 3.4). As the glacier outlines are derived at the end of the ablation season, the climate parameters are aggregated over the preceding winter and the related summer months. The significance of the correlation coefficient with a confidence level of 90% was estimated by the help of the student distribution according to Obilor and Amadi (2018).

Regression of Glaciated Area
The mean change of the glaciated area per elevation zone is estimated via a least-square robust adjustment using a linear model. Concerning the decadal trends observed in the meteorological data (cf. Fig. 2), the trend is analyzed separately for each decade after 1990, as a sufficient number of images have been available since then. The significance of the observed trends is estimated by the help of the Mann-Kendall test with a confidence level of 90%. Additionally, the total trend over the whole time frame is taken into account: Fig. 4 depicts a square gradient, and Tables 3, 4, 5 lists the derived linear gradients in the last column.

Prediction of Glaciated Area
From the trends observed over the whole time period, future behaviour can be estimated with a measurable uncertainty.
The linear model is utilized to provide a prediction for each elevation zone under study (cf. Fig. 7). To consider other mathematical models as well, the linear model from Fig. 7 and the quadratic model from Fig. 2 are completed with an The ArcticDEM is used to display the glacier boundaries from Fig. 3 according to their elevation exponential model. These models are applied to the total glaciated surface area exclusively for stability reasons. The respective degree of determinism-similar to the correlation coefficient introduced above-and the standard deviation of the prediction are reported to describe the goodness of fit. Finally, the glaciated area is predicted for the coming decades with the upper and lower limit of the 90% confidence interval in Table 6.

Equilibrium Line Altitude
The ELA is estimated via the snowline on images recorded at the end of the ablation season along with selected profiles. The NDSI is used to derive the glaciated areas for 1973 (all three glaciers), 2016 (Langjökull), 2017 (Vatnajökull) and 2018 (Hofsjökull). To derive the accurate boundary of the accumulation (snow) and ablation (ice) zones, a maximum likelihood classification (Bronge and Bronge 1999) was carried out. The classes snow, firn/ice, and tephra cover can be differentiated (Dozier 1989) due to their spectral properties (cf. Table 1). Since the tephra-covered areas lie exclusively within the ablation zone, those areas can be merged with the firn/ice class. The equilibrium line altitudes are then calculated by the intersection of the end-of-summer snowline with the ArcticDEM along the profiles drawn in Fig. 1.

Results
The meteorological data are illustrated in Fig. 2. It highlights the mean summer temperature and total winter precipitation in Stykkishólmur since 1949. A slight increase in the mean summer air temperatures can be noted. The graph highlights certain peaks and lows, whereby the lows are noted at the beginning of the study period, i.e., 1973 with 8.45 °C mean summer air temperature, and the peaks in the more recent years, i.e., 2010 with 11.08 °C, 2014 with 10.63 °C or 2016 with 10.35 °C. Regarding the total winter precipitation, a slight increase from the 1940s to now can be detected as well. The early 1990s-especially 1991 and 1993-highlight the years under study with the highest observed winter precipitation. The winter precipitation is declining with the lowest precipitation in winter 1998/99. An increase in temperatures is contrasted with a slightly wetter winter in the following year. This, in particular, can be observed in

Glaciated Area
In Figs. 3, 4, the area changes of the Langjökull glacier from 1941 until 2018 are depicted. The glacier is evidently changing steadily, meaning it is growing and shrinking over time. From a long-term view, however, there is a decline in the glacier surface area. The Hofsjökull glacier has its maximum extent on its north-eastern side around 1993. On its north and northwestern flank, the glacier experienced its largest expansion in 1986. A general retreat of the glacier from the turn of the millennium on can be detected. Vatnajökull shows a strong decline in the first decade of the new millennium after a decade of relative stability in the 1990s. Since 2010, the glaciated surface is increasing again. All three plots indicate that the increasing mean summer temperature (cf. Fig. 2) is related to a stronger negative gradient of the glaciated area. Heavy short-term fluctuations in the mean summer temperature evoke a direct response of the glaciated surface (cf. Figs. 2, 4) from 2008 on. However, constant temperatures do not automatically lead to a steady glaciated area (cf. Figs.2, 4 from 1991(cf. Figs.2, 4 from until 2002. As indicated by Figs. 2, 4, a decrease of the glaciated surface area can also be observed in years of lower temperatures. The diagrams in Fig. 4 only depict the total surface of the respective glaciers. The parameter elevation, which is known to have a significant influence (Raup et al. 2007), is not yet considered.

Glaciated Area by Elevation
The multivariate histograms in Fig. 5

Correlation with Climate Parameters
To transfer the observations by time to the climate parameters, the glaciated area by elevation zone is correlated with the measurements at Stykkishólmur station. Fig. 6 illustrates the derived Pearson correlation coefficients: blue for negative and red for positive correlation, grey indicates a nonsignificant correlation. The elevation zones form columns, the selected climate parameters are displayed in rows. The decrease with time explored in the preceding section is confirmed for lower elevation zones on all three glaciers, cf. blue boxes in the first row of Fig. 6. This correlates almost equally strong with the mean winter and even more with the mean summer temperature. Only Vatnajökull shows a different behaviour: exclusively the elevation zones between 700 m and 1200 m show a negative correlation with the mean summer temperature. The same elevation zones are subject to a decrease with increasing summer precipitation. This interrelation can also be observed on Langjökull and Hofsjökull. They also show a slight negative, but significant, correlation with the air pressure in summer. Another negative influence is registered for sunshine hours. Cloud cover, in contrast, has a positive effect both in winter and in summer. The most surprising correlation is the dependence of the glaciated area on the mean winter wind speed at both Hofsjökull and Langjökull for all altitude zones below 1400 m. Concerning Vatnajökull, this behaviour is only significant between 800 m and 1000 m. Vatnajökull surprises with a positive correlation to the mean summer temperature, summer precipitation, and air pressure in summer in-between 1500 m and 1900 m in elevation. Other positive trends are recorded for winter cloud cover on the lowest parts of the glacier and summer wind speed for the elevation zone of 1800 m. In comparison, the correlations observed at the elevation zones covered by all three glaciers (i.e., 700-1500 m) are very similar.

Regression of Glaciated Area
The decadal trends of the glaciated surface area, according to the trends estimated from the summer temperature and the winter precipitation, are shown in Tables 3, 4, 5. For Hofsjökull in Table 3, the elevation zones above 1500 m do not show any significant trend. Regarding the latest time span 2010-2018, it is even assured that there is no observable trend. For the whole period 1950-2018, the elevation zones above 1200 m record only negligible trends. By contrast, the zones from 900 m to 1200 m are characterized by high variability with an increase from 1950 to 1990 and from 2010 to 2018 and a decrease in all other time spans. Lower elevation zones show a steady decrease with its maximum in the first decade of the new century with a loss of 122 ha per year, which equals − 0.57%. In sum, a significant negative trend can be stated for all elevation zones below 1200 m that reaches almost − 1% per year.
The glaciated area of Langjökull has a similar pattern. Table 4 lists only minor to negligible trends ( < 6 ha ) for elevation zones above 1200 m except for the period between 1990-2000, when a decrease of 12 ha was recorded. In this Fig. 6 Correlation of the glaciated area in different elevation zones and selected climate parameters. Only significant correlations according to a 90% confidence intervall are depicted in colours. Non-sig-nificant correlations are grey. Unavailable information due to smaller vertical extension of the respective ice cap is symbolized by the white background period of time the elevation zones between 800 m and 1100 m are subject to an extreme decrease of more than 1 km 2 per year. The subsequent decade 2000-2010 shows a slight increase for the elevation zones 1000 m to 1200 m, but overall, the steady decrease is dominant. Elevation zones below 900 m show a consistent decrease, whereby the lowest value is − 0.74 km 2 per year from 2000 on.
With view to Vatnajökull in Table 5, the large vertical extent of the ice cap from almost sea level to 2200 m is noticeable. The change values are higher because of the larger area covered by Vatnajökull. Therefore, the relative changes are also listed. In contrast to the preceding two examples, Vatnajökull is characterized by a high variability over almost all elevation zones. Only zones below 300 m show a steady decrease reaching 4.5% in the decade 1990-2000 and even 5.8% in the first decade of the new century. Regarding the whole study period, minor or nonsignificant positive trends are only observed for elevation zones between 600 and 700 m and above 2000 m. All other zones show a decrease of at least 1‰ per year, increasing to 1.2% towards the lowest parts of the glaciers. Figure 7 shows the applied linear regression with the prediction per elevation zone. Obviously, the lower parts of each glacier will disappear or reduce extremely in the coming decades. Since short-term fluctuations-as visible for the time span under study-are not considered in the prediction model, the predicted decrease is linear until the respective elevation zone is completely ice-free.

Prediction of Glaciated Area
To verify the suitability of the regression model, other models are also taken into account. The parameters derived from linear, quadratic and exponential models are used to predict the change in total glacier area up to the year 2100 (cf. Table 6). The best-fitting model for Hofsjökull is the exponential model with a determinism coefficient of almost 80% and a standard deviation of 30 km 2 . It only slightly excels the linear model displayed in Fig. 7. With respect to Langjökull, the quadratic model delivers the most reliable predication with a certainty of more than 80% and only 22 km 2 standard deviation. But, Gradients that did not pass the Mann-Kendall test with a confidence level of 90% are greyed out Table 4 Adjusted increase or decrease in hectares (left number) or per thousand (right number) per year for different levels and over varying time spans for Langjökull Gradients that did not pass the Mann-Kendall test with a confidence level of 90% are greyed out it predicts that the glacier will have disappeared until 2100. In the optimal case, it only reduces to a remaining size of 254 km 2 according to this model. Although the quadratic model delivers the best results for Vatnajökull, the goodness of fit is similar over all three models.
Accepting that the median of the three predictions listed in Table 6 is the most likely case, we can state a size of only 338 km 2 for Hofsjökull in the year 2100, i.e., only 40% of its current extent. Langjökull will reduce to 700 km 2 , which still equals 80% of the present extent. Vatnajökull, although highly variable in time, could recede to Gradients that did not pass the Mann-Kendall test with a confidence level of 90% are greyed out Fig. 7 Temporal change of glaciated areas per elevation zone from 1949 until 2018 (saturated colours) including a linear prediction (leach-out colours) for three further decades based on the whole time frame under study 6279 km 2 , which is also about 80% of its current size. The most negative model for the glaciers predicts a complete disappearance of Hofsjökull and Langjökull and a reduction of Vatnajökull to 6% compared to its extent in 2017.

Equilibrium Line Altitude
The EL for 1973 is shown in dark blue, while the most recent EL (Langökull 2016;Hofsjökull 2018;Vatnajökull 2017) is shown in turquoise. The ELA measured along selected outlet glaciers are shown in Table 7 (i.e., "L1": glacier Langjökull, profile 1). The accumulation area decreased significantly on all three glaciers. Thus, the mean increase over this period is more than 200 m on all three glaciers studied and the minimal increase varies from 38 m to 111 m. The maximum increase starts at 360 m and reaches 592 m in the case of Vatnajökull.

Glaciated Area
To derive the glacier outlines, the definition of a threshold value is required above which a pixel is classified as "snow" or "ice", respectively. This step must be done carefully to prevent an over-or underestimation of the snowand ice-covered surface. The threshold value differs across the selection of sensors. The lower threshold values used for extracting the glacier bodies range between 0.35 and 0.51 for the Landsat 5, 7, 8 and Sentinel-2 images. The values are close to those reported in the literature for the NDSI, namely between 0.25 and 0.45 (Hall et al. 1995). Only the threshold range of aNDSI derived from the Landsat 1 images differs slightly [0.015 − 0.370] as expected (cf. Sect. 4.1). One difficulty here is to set a suitable acceptance range instead of only a lower threshold value, which completely includes the Skeiðarárjökull outlet glacier in the South of Vatnajökull. This glacier is partially covered with tephra layers (pyroclastic sediment) since this glacier is located near an active volcano. The glacier sites overlaid with tephra layers are not well recognized by spectral indices. Due to their spectral properties, tephra can be mistaken for rock or soil. They are therefore not included in the glacier mask. There are different ways to derive such covered glaciers. One approach follows the method of combining band ratios (Landsat 5 TM bands 3, 4 and 5) and the use of a slope map (Paul et al. 2004). Since the Skeiðarárjökull glacier tongue is very shallow, the derivation by a combination of a Digital Elevation Model (DEM) and satellite data, is not possible. Therefore, a sufficiently large threshold is chosen that includes both the glaciers as well as the Skeiðarárjökull outlet glacier. Figure 4 shows the changes of the glaciated surface area over time for the three glaciers studied with decadal linear trends and a quadratic regression over the entire time span. Non-significant trends are greyed out. Langjökull, Hofsjökull, and Vatnajökull cover-according to historical maps from the 1940s-an area of 985 km 2 , 952 km 2 and 8280 km 2 respectively. While both Langjökull and Hofsjökull show a retreat in the years 1941-1973, Vatnajökull even seems to grow slightly. The average summer temperatures show a negative trend in the same period. While winter precipitation indicates quite high amounts, especially in the period from 1949 to 1960, snowfall amounts decrease in the following decade until 1972 (cf. Fig. 2).
In contrast to Langjökull and Hofsjökull, where a slight glacier growth is observed between 1973 and 1986, Vatnajökull stagnates for almost ten years. However, in the following years, a general glacier retreat can also be observed at Vatnajökull. Although total winter precipitation and summer temperatures continue to increase linearly, strong negative outliers are observed, especially with regard to summer temperatures. In the following years, the positive temperature trend continues. Langjökull and Hofsjökull simultaneously lost glaciated areas of a significant amount until 1991. Particularly for Hofsjökull, a glacier area loss of almost 100 km 2 is recorded. Both Hofsjökull and Langjökull increase in size during the period 1991-1993 (cf. Fig 4). In parallel, significantly lower summer temperatures, with a decrease of about 1 °C compared to previous years, and increased snowfall were recorded at the Stykkishólmur Weather Station. Over the next decades, until 2018, all three glaciers show an overall negative trend in terms of glaciated area. In those years, a clear positive trend can also be seen with regard to the average summer temperatures. With regard to winter precipitation, increased snow amounts are indicated, but significant outliers are also evident (cf. Fig 2). In 2000 the glaciated area of the ice cap of Langjökull is 926 km 2 , the area of Hofsjökull is 894 km 2 and Vatnajökull covers an area of 8108 km 2 . Björnsson and Pálsson (2008) determine for the glaciated area of Langjökull in the 2000s 900 km 2 , for the area of Hofsjökull 890 km 2 and Vatnajökull 8100 km 2 , which equals a deviation from our estimation of less than 3% in the case of Langjökull and even less than 1% in the case of Hofsjökull and Vatnajökull. The negative trend in glacier area continues, especially at Langjökull and Hofsjökull, but also at Vatnajökull until 2018. At Vatnajökull, however, there is also a recognizable and significant increase in the glaciated area from 2014 to 2015. Pálsson et al. (2015), studied the variation in mass balance at Vatnajökull 2014-2015, noting  1973-2016 1973-2018 1973-2017 that the 2014-2015 glacier year is the first in two decades with positive mass balance for the glacier.

Glaciated Area by Elevation
The area change at both Hofsjökull and Langjökull occurs exclusively at elevations below 1300 m (cf. Fig. 5). Lower elevations (< 1300 m at Hofsjökull and < 1100 m at Langjökull) show a high variation over time, while higher elevations show no variation at all. Though the slight increase (in green) is visible, the trend indicates a steady decrease, cf. the chronological order of blue, green, and orange at the left edge in Fig. 5. The distribution of glaciated areas according to the elevation zone is quite different for Vatnajökull in comparison to the other glaciers under study. Fig. 5 shows a vast decrease at elevations around 1400 m and a significant increase at elevations around 1600 m around the year 1990. Hence, all elevation zones are subject to high variations. Only a few parts of the glacier exceeding 2000 m in elevation remain relatively constant. The main reason for the observed variation at higher elevation is the volcanic activity causing ash and tephra cover, so that the underlying ice is no longer visible for remote sensors.

Correlation with Climate Parameters
The correlation analysis of the glaciated area in different elevation zones and the respective climate parameters such as mean temperature, total precipitation, air pressure, cloud cover, and wind speed for summer (June-September) and winter (October-April) based on the Stykkishólmur weather station is provided in Fig. 6. The analysis indicates especially on Langjökull and Hofsjökull but also on Vatnajökull an exceptionally high correlation of glacier growth with an increased winter wind speed. Colder summer as well as winter temperatures, more clouds in summer and thus fewer hours of sunshine also correspond with increased glacier growth. It should be noted that the Stykkishólmur measuring station is neither close to the three observed glaciers (cf. Fig. 1), nor is the difference in altitude taken into account. Therefore, it is quite possible that climate parameters might differ locally from the ones measured at Stykkishólmur. Nevertheless, the meteorological measurements used here are the best available ones for the region of interest (cf. Sect. 3).

Regression of Glaciated Area
The decadal trends of the glaciated area by elevation confirm the trends estimated from summer temperature and winter precipitation (Tables 3,4,5). They indicate constant extension at higher altitudes, but continuous shrinkage at lower elevations. Vatnajökull is the only glacier characterized by several positive trends in higher regions, which can be referred to volcanic activities in the first decade of the new century. Nevertheless, the total glaciated surface area is also shrinking (cf. Fig. 7).

Prediction of Glaciated Area
The estimated long-term development is depicted in Fig. 7. For this purpose, three models, linear, quadratic, and exponential, are being tested with expectation value and upper and lower 90% limits. The quadratic model delivers the most reliable prediction with a certainty of more than 80% and only 22 km 2 standard deviation for Langjökull. The largest extension is reached in 1993 (cf. Fig. 4). The higher parts of the glacier are quite stable except for two years around 2010. Glaciated areas below 500 m have already almost disappeared. The areas between 500 and 1000 m are continuously shrinking.
For Hofsjökull, the exponential model with a certainty of almost 80% and a standard deviation of 30 km is considered the most appropriate. The expansion at higher altitudes is constant, but a continuous shrinkage at lower altitudes has been demonstrated. Two local peaks of maximum expansion can be observed in the 1980s and 1990s, and a local minimum in 2010, resulting from a decrease in all altitude zones below 1500 m. This is remarkable because most deviations are only observed in lower regions so far.
All tested models fit equally well for Vatnajökull. Numerous local extrema can be observed since the 1980s. These variations can be referred to several volcanic eruptions in recent years. At second glance, an interesting detail is additionally depicted: The glaciated areas below 1000 m are shrinking and the areas above 1500 m do not change, i.e. only the areas between 1000 m and 1500 m, cf. the turquoise polygon in Fig. 7, are expanding. Unfortunately, this detail cannot be explained in a satisfactory manner based on the present data set and has to be followed in future studies.

Equilibrium Line Altitude
Additional information on the ELA variations between 1973 and 2018 are gathered using selected profiles as seen in Fig.1. The ELA of the Langjökull migrates significantly higher between 1973 and 2016. The largest retreat on Langjökull outlet-glaciers is noted along the profiles "L1", "L4", "L9", and "L10" (Fig.1 and Table 7). These have withdrawn by more than 30% in comparison to their previous location in 1973. An extraordinary fluctuation is noted on profile "L9", which moved close to 400 m higher. This corresponds to a retreat of rounded 50% from 2016 to 1973. Hofsjökull's ELA is observed between 1973 and 2018. On its outlet glaciers, located at "H9" and "H11", the ELA has withdrawn the furthest ( Fig. 1 and Table 7). Profile "H11", which corresponds to Hofsjökull's northwestern flank, has retreated more than 360 m during the investigation period (Table 7).
On the Vatnajökull glacier, the ELA is evaluated between 1973 and 2017. Particularly striking is the ELA along profile "V4" and even more along "V5" (Fig. 1 and Table 7). At the latter location, "V5", the ELA has increased by 592 m. On average, the ELA at Vatnajökull retreats by 205 m (Table 7). Björnsson and Pálsson report an ELA fluctuation from 200 m-300 m during the study period 1992-2007 (Björnsson and Pálsson 2008). Guðmundsson et al. (2017) studied changes of Breiðamerkurjökull, south-eastern Vatnajökull, from its late nineteenth-century maximum to the present. They reported an estimated ELA in the 1890s of between 850 m-950 m. The location of Breiðamerkurjökull corresponds to profile "V6" in this study, where an ELA of 934 m in 1973 was assessed.

Conclusion
Based on Landsat 1, 5, 7, 8 and Sentinel-2 imagery from 1973 to 2018 as well as data obtained from historical maps from the 1940s, we compared glacier termini fluctuations of the Langjökull, Hofsjökull and Vatnajökull glaciers in Iceland to registered climate parameters. All three ice caps retreated during the investigation period with a particularly sharp decline since the turn of the millennium. Since the 1940s, the Langjökull glacier has decreased by more than 100 km 2 , the Hofsjökull over 115 km 2 and the Vatnajökull over 430 km 2 . Area changes at Langjökull glacier exclusively occur at elevations below 1300 m . Lower elevations are characterized by a high variation over time, but indicating a steady decrease, whereas glaciated areas below 500 m have already almost disappeared. This steady decline can also be seen at the Hofsjökull glacier. A retreat is generally detectable at the Vatnajökull, but the distribution of the glaciated areas according to the elevation zone additionally depends on local impacts such as volcanic activities. The analysis of climate data for Iceland indicates that ice-frontal advance can be associated with an increased wind speed in winter, colder summers as well as lower winter temperatures, more clouds in summer and therefore fewer sunshine hours. To estimate a long-term development, three models (linear, quadratic and exponential) are tested with the expectation (mean) value as well as upper and lower boundaries. At Langjökull, glaciated areas below 500 m have already almost disappeared, and at Hofsjökull, the extension at higher elevations is constant, but a continuous shrinking at lower elevations is proven. With respect to the coming years, the steady retreat is most likely (cf. Fig. 7). ELA variations on several outlet-glaciers of Langjökull, Hofsjökull and Vatnajökull were assessed, by obtaining imagery at the end of mass-balance years 1973 and 2016/17/18. The ELA migrates significantly higher on all surveyed outlet-glaciers.
This study highlights not only the need to place observations of contemporary glacier changes in a broader, longerterm context, but also the great use of remote sensing data for continuous monitoring of ice caps and glaciers. The results could be further improved by the use of satellite data with higher spatial and temporal resolution. On-site measurements could possibly monitor the tephra-covered glacier tongues with the highest precision to include small and overshadowed glacier parts. Optical remote sensing is limited in terms of cloud cover and solar illumination, which both impede the use of optical satellite images over Iceland. Radar remote sensing would have the possibility to penetrate clouds and to operate independently from sunlight, but, the sensitivity to dry snow cover is rather limited. Therefore, a synergy of optical sensors and radar systems is needed for a comprehensive glacier observation as stated in literature (Nagler et al. 2016). Though the results derived from 2D image data coincide very well with the locally evaluated mass balances, the temporal change of height still is one key information to be included in future studies to assess the mass balance of glaciers with high spatial and temporal resolution.