Glacier Changes in Iceland From ∼1890 to 2019

The volume of glaciers in Iceland (∼3,400 k m 3 in 2019) corresponds to about 9 mm of potential global sea level rise. In this study, observations from 98.7% of glacier covered areas in Iceland (in 2019) are used to construct a record of mass change of Icelandic glaciers since the end of the 19th century i.e. the end of the Little Ice Age (LIA) in Iceland. Glaciological (in situ) mass-balance measurements have been conducted on Vatnajökull, Langjökull, and Hofsjökull since the glaciological years 1991/92, 1996/97, and 1987/88, respectively. Geodetic mass balance for multiple glaciers and many periods has been estimated from reconstructed surface maps, published maps, aerial photographs, declassified spy satellite images, modern satellite stereo imagery, and airborne lidar. To estimate the maximum glacier volume at the end of the LIA, a volume–area scaling method is used based on the observed area and volume from the three largest ice caps (over 90% of total ice mass) at 5–7 different times each, in total 19 points. The combined record shows a total mass change of −540 ± 130 Gt (−4.2 ± 1.0 Gt a − 1 on average) during the study period (1890/91 to 2018/19). This mass loss corresponds to 1.50 ± 0.36 mm sea level equivalent or 16 ± 4% of mass stored in Icelandic glaciers around 1890. Almost half of the total mass change occurred in 1994/95 to 2018/19, or −240 ± 20 Gt (−9.6 ± 0.8 Gt a − 1 on average), with most rapid loss in 1994/95 to 2009/10 (mass change rate −11.6 ± 0.8 Gt a − 1 ). During the relatively warm period 1930/31–1949/50, mass loss rates were probably close to those observed since 1994, and in the colder period 1980/81–1993/94, the glaciers gained mass at a rate of 1.5 ± 1.0 Gt a − 1 . For other periods of this study, the glaciers were either close to equilibrium or experienced mild loss rates. For the periods of AR6 IPCC, the mass change rates are −3.1 ± 1.1 Gt a − 1 for 1900/01–1989/90, −4.3 ± 1.0 Gt a − 1 for 1970/71–2017/18, −8.3 ± 0.8 Gt a − 1 for 1992/93–2017/18, and −7.6 ± 0.8 Gt a − 1 for 2005/06–2017/18.


INTRODUCTION
Glaciers in most areas of the world are losing mass as global temperatures rise in response to increased greenhouse gas concentrations in the atmosphere (e.g., Vaughan et al., 2013;Hock et al., 2019;Meredith et al., 2019;Zemp et al., 2019). In situ mass-balance observations are sparse (e.g., Zemp et al., 2020a), but with the aid of satellite and other remote-sensing data, an increasingly clear picture of glacier mass loss around the world has been appearing (e.g., Brun et al., 2017;Wouters et al., 2019;Morris et al., 2020;Shean et al., 2020). Glacier mass loss is a global phenomenon, and the rates in the early 21st century are unprecedented for the observed period (Zemp et al., 2015). Reconstructions of glacier mass-change rates for the 20th century and the first decade of the 21st century show substantial temporal and spatial variations, but a global mass loss trend became clear toward the end of the 20th century (Leclercq et al., 2011;Marzeion et al., 2015;Marzeion et al., 2012).
Iceland is located in a region of maritime climate in the middle of the North Atlantic Ocean with relatively cool summers, mild winters, and high precipitation. Glaciers in Iceland are all temperate and cover about 10% of the area of the country (Björnsson and Pálsson, 2008), with the largest ice cap Vatnajökull (∼7,700 km 2 , ∼2,870 km 3 , in the year 2019) located near the southeast coast, two smaller ice caps Langjökull (∼835 km 2 , ∼171 km 3 , in the year 2019) and Hofsjökull (∼810 km 2 , ∼170 km 3 , in the year 2019) in the central highlands [area estimates are from Hannesdóttir et al. (2020) and volumes are calculated in this study], and Mýrdalsjökull [∼598 km 2 , ∼140 km 3 , in the year 1991 (Björnsson et al., 2000)], near the southern coast. Seven additional glaciers in Iceland are larger than 10 km 2 , and there are presently around 250 smaller glaciers, many of them in the central north highlands (Tröllaskagi (Trö) in Figure 1). The inventory of Icelandic glaciers made around the year 2000 includes about 300 glaciers (Sigurðsson and Williams, 2008). An update of this inventory in 2017 showed that some tens of those had disappeared or were categorized as dead ice . Several glaciers had, during their retreat, split into two or more separate glaciers or ice patches. The area loss since the end of the Little Ice Age (LIA) is ∼2,200 km 2 and ∼750 km 2 since the year 2000, or about 40 km 2 (or 0.4%) per year (Hannesdóttir et al., 2020).
Glaciers in Iceland have received much attention through the centuries due to their proximity to inhabited regions ( Figure 1). The retreat of many glacier tongues was noticed in the early 20th century and in 1930 a country-wide voluntary monitoring program was initiated (Eyþórsson, 1963;Sigurðsson, 2005). This program was later continued by the Icelandic Glaciological Society (founded in 1950) and continues to this day (Björnsson, 2017;Hannesdóttir et al., 2020). Thirty to forty terminus positions are measured annually and the observations are posted on the website spordakost.jorfi.is. The surface and bedrock topographies of the largest ice caps have been measured in radio-echo sounding campaigns carried out since 1977. Radioecho sounding on the temperate ice caps in Iceland required a much longer electromagnetic wavelength than had been used on the cold polar ice caps (Björnsson and Pálsson, 2020). The pioneering measurements resulted in a good knowledge of the subglacial topography of the ice caps and their total volume FIGURE 1 | Map of Iceland showing the glaciers considered in this study. The pie chart in the bottom right corner shows the relative sizes of the glaciers in 2017 (Hannesdóttir et al., 2020) in percentage; Vatnajökull (V), Langjökull (L), Hofsjökull (H), and the smaller glaciers grouped as O m in pie chart: Drangajökull, Snaefellsjökull, Eyjafjallajökull, Tindfjallajökull, Torfajökull, Mýrdalsjökull, Þrándarjökull, Hofsjökull eystri, Tungnafellsjökull, Hrútfellsjökull, Eiríksjökull and Barkárdals-and Tungnahryggsjökull, considered as one glacier unit on Tröllaskagi. Unmeasured glaciers (O u ) corresponded to ∼1.3% of the total area in 2019 (Hannesdóttir et al., 2020). The town Stykkishólmur is shown with a purple dot, from which a temperature record exists since the middle of the 19th century.
Mass loss from Icelandic glaciers in the years 1994/5-2009/10 was reported as 9.5 ± 1.5 Gt a −1 with a large interannual variability . In this article, we present an updated and extended record for the hydrological years 1890/ 91 to 2018/19. The main objective was to derive an estimate for the mass-change history of Icelandic glaciers for the 20th century and the first two decades of the 21st century. All recent studies and measurements with various methods are compiled and combined; a new estimate of the ∼1890 volumes of the three largest ice caps is made using new data on their area at that time (Hannesdóttir et al., 2020) and volume-area scaling (Bahr et al., 2015), and the recently established non-surface mass balance of Icelandic glaciers, that takes into account energy dissipation caused by the flow of water and ice, geothermal melting, volcanic eruptions, and calving , is included to improve the estimate of the mass-balance history of all glaciers in Iceland since the end of the LIA.

METHODS AND DATA
The time steps in this study correspond to the glaciological year  from autumn to autumn, using the floatingdate mass-balance system (Østrem and Brugman, 1991;Björnsson et al., 2002;Cogley et al., 2011), that is, the end of the summer melt season marks the start of a new glaciological year. To construct the mass-balance history of Icelandic glaciers back to the year 1890/91 [assumed here to be the end of the LIA in Iceland (Thorarinsson, 1943;Sigurðsson, 2005)], we compile and combine different data sets and methods. These include the following: (1) Record of the annual surface mass balance obtained with the glaciological method for the three largest ice caps, Vatnajökull (glaciological years 1991/92 to 2018/19), Hofsjökull (1987/88 to 2018/19), and Langjökull (1996/97 to 2018   and is forced with a ERA-Interim downscaling using the HARMONIE-AROME model at 2.5 km resolution (Schmidt et al., 2019) (see shaded area in Figure 3A). Data from automatic weather stations and glaciological surface mass balance, and runoff measurements were used to constrain the model (Schmidt et al., 2018). (3) The non-surface mass-balance estimates from Jóhannesson et al. (2020) are taken into account and added to the glaciological and modeled surface mass balance listed above. (4) Zero annual mass balance for Vatnajökull for the period 1970/71 to 1979/80 and for Hofsjökull from 1970/71 to 1986/ 87 is assumed (green boxes in Figures 3A,B). This is supported by the small area changes during these periods (Hannesdóttir et al., 2020), the relatively cold weather conditions ( Figure 3E), the observations of limited terminus length changes (spordakost.jorfi.is), and the near-zero geodetic mass balances of several of the smaller glaciers in these periods ( Figure 3D). (5) Geodetic mass-balance records for Langjökull (Pálsson et al., 2012), from 1937/38 to 1996/97 (red lines with uncertainties in Figure 3B) and 12 smaller glaciers ( Figure 3D) from 1945/ 46 to 2016/17 ) that cover 8.3% of the glacier area in Iceland. The remaining 1.3% of glacier area that is not observed is assumed to have the same mass loss rate as the measured small glaciers. The geodetic results were also used to scale  and validate (Pálsson et al., 2012) the glaciological measurements. (6) Estimates of the volumes of Vatnajökull and Hofsjökull in 1890 and 1945 as well as the volume of Langjökull in 1890 based on the volume-area scaling (Bahr et al., 1997;Bahr et al., 2015) ( Figure 4). The glacier areas are derived from Hannesdóttir et al. (2020); we describe the method below and the resulting mass change rates calculated from the area and volume changes are shown with purple boxes in Figures 3A Belart et al. (2020). After 1994/95, an estimate of the annual variability of the mass change for other glaciers than the three largest is included by calculating the net mass change for each glaciological year, using the corresponding F value for each period. The net mass change during these periods, which is obtained with the geodetic method, is not altered by this.
The above records were combined to calculate the mass balance of all Icelandic glaciers from 1890/91 to 2018/19. The cumulative mass change of all the glaciers for each period was computed as the sum of the total mass balance of the all glaciers (the product of the specific mass balance and the time-dependent glacier area). Below we discuss i) how the surface mass balance from the glaciological method is combined with the non-surface mass-balance estimates, ii) the application of the volume-area scaling to estimate past volumes of the three largest ice caps, and iii) the uncertainty of the obtained total mass balance of the Icelandic glaciers.

Combining the Surface and Non-Surface Mass-Balance Records
The surface mass balance from the glaciological method is obtained by measuring the snow water equivalent (w.e.) thickness of the snowpack above last summer's surface in spring, and the summer ablation at the end of the summer melt season at ∼60 survey sites on Vatnajökull and ∼25 survey sites on Hofsjökull and Langjökull. These measurements, together with the HIRHAM5 snowpack model simulations, yield a detailed record of the winter, summer, and annual surface mass balance during the last three to four decades for these ice caps (Figure 2), constituting ∼90% of the glacier area in Iceland. These records show rather large annual variability and the exceptional effect of the 2010 eruption in Eyjafjallajökull, which doubled the summer melt due to the deposition of volcanic tephra on the surface of the ice caps .
The non-surface mass balance is known to be a non-negligible component of the mass balance of glaciers in Iceland (e.g., Björnsson et al., 2013) but has so far not been included in mass-balance estimates. The previously published massbalance record for the Icelandic glaciers has now been revised by including this component. The non-surface melting of glaciers in Iceland for the period 1995-2019 estimated by Jóhannesson et al. (2020) includes geothermal melting, energy dissipation caused by the flow of water and ice, volcanic eruptions, and calving.
For Hofsjökull and Langjökull, the non-surface mass balance due to volcanic eruptions and calving is negligible, but not for Vatnajökull. Rather than including the ice melt (∼3.7 Gt) due to the Gjálp eruption in October 1996 (Guðmundsson et al., 2004), in the average rate over the period (1995/95-2018/19), it is subtracted and then added to the glaciological year 1996/97. This reduces the average mass balance due to geothermal melting and volcanic eruptions on Vatnajökull presented by Jóhannesson et al. (2020) from −0.075 m w.e. a −1 to −0.055 m w. e. a −1 . The mass loss due to energy dissipation in Vatnajökull, caused by the flow of water and ice as estimated by Jóhannesson et al. (2020), is assumed constant for the whole period: −0.085 m w.e. a −1 .
Mass loss due to calving is negligible for all outlets except for Breiðamerkurjökull on the south side of Vatnajökull, which calves into the terminal lake Jökulsárlón. We apply a variable calving rate as described by Jóhannesson et al. (2020) for the years 1996/97-2016/17, resulting in an average of −0.056 m w.e. a −1 .
For the periods without observations, we assume a calving rate that changes linearly from 0 in 1950/51 [start of significant calving (Björnsson et al., 2001)] to −0.029 m w.e. a −1 in 1996/ 97 and use −0.067 m w.e. a −1 (value in 2016/17) for the last two years of the record.
The average non-surface mass-balance values for Langjökull and Hofsjökull, where calving is insignificant and no eruption has taken place during our study period, are −0.055 m w.e. a −1 and −0.067 m w.e. a −1 , respectively. In Figures 3A,B,C, the nonsurface mass-balance estimates have been added to the annual surface mass-balance records.
Comparison of the glaciological surface mass-balance record of Hofsjökull with results from geodetic mass balance, derived by differencing digital elevation models (DEMs), revealed a bias between the two data sets. This bias has been corrected in the surface mass-balance record of Hofsjökull Thorsteinsson et al., 2017) shown in Figure 2, taking into account the contribution of the non-surface mass balance.

Volume-Area Scaling
The volume of a glacier can be obtained by integrating the ice thickness, calculated as the difference between surface and FIGURE 2 | Mean specific surface mass-balance records of the three largest glaciers in Iceland obtained with glaciological method (in situ surveys). The mass-balance measurements have been conducted at ∼60, ∼25, and ∼25 locations since 1991/92, 1996/97, and 1987 (Korona et al., 2009;Pálsson et al., 2012). The areal extent at several times for all the glaciers has been estimated based on a data set of glacier outlines (Pálsson et al., 2012;Jóhannesson et al., 2013;Hannesdóttir et al., 2020). Ice-volume estimates at other times can be calculated by multiplying the annual specific mass balance ( Figure 3) of each glacier by the corresponding glacier area, linearly interpolated with time between dates of area observations, converting the annual mass change into ice volume [assuming the conversion factor 0.85 (Huss, 2013); note that mass-balance records previously published that used conversion factor 0.9 (Pálsson et al., 2012;Jóhannesson et al., 2013) have been adjusted accordingly (Thorsteinsson et al., 2017)] and integrating the volume change relative to the date of the surface DEMs listed above. This calculation gives the area and volume for several times back to 1970 for Hofsjökull and Vatnajökull and to 1937 for Langjökull. The results for the three ice caps for dates of observed areas and corresponding glacier volumes are shown with filled circles on the volume-area scatter plot in Figure 4. Volume-area scaling, as described by Bahr et al. (1997), Bahr et al. (2015), gives an empirical relation between glacier volume (V) and area (A) of the form: This relationship was already applied to estimate the ∼1890 volume of Langjökull by Pálsson et al. (2012). Here, we assume that the relationship between area and volume is the same for the three largest glaciers in Iceland, Vatnajökull, Langjökull, and Hofsjökull, which is supported by the linear relationship shown in Figure 4. The available area and volume data for all of them are therefore combined to estimate the parameters for the volume-area scaling equation. We subtract the area of debris-covered moraines and dead-ice fields with little surface-elevation change in the GLIMS glacier outlines (Hannesdóttir et al., 2020) (in total 27 km 2 in the year 2019), Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 5 in order to base our volume-area relationship on the area of the active glacier where substantial ice-thickness changes take place (these debris-covered areas are specified as polygons in our GLIMS data set). Using the least-squares fitting method, we found c 0.0379 and c 1.25; the value of the exponent c is the same as obtained independently by Bahr et al. (1997) for ice caps.
Volume-area scaling has been widely used to estimate the volume of glaciers with an unknown subglacial topography/ice thickness (e.g., Radić and Hock, 2010), and changes in ice volume associated with variations in glacier extent when multi-temporal DEMs of the glacier surface are not available (e.g., Pálsson et al., 2012). The method is found to be uncertain by tens of per cent, up to even a factor 2-3, for estimating the volume of ice caps with an unknown subglacial topography (Gärtner-Roer et al., 2014), and methods that include information on glacier mass balance and glacier-surface geometry using ice-flow dynamics (Huss and Farinotti, 2012;Farinotti et al., 2019) are preferred for this purpose. Using volume-area scaling to estimate changes in the volume of glaciers with a well-known subglacial topography, from variations in glacier area over decadal time spans, may be expected to be more accurate because this mainly relies on the assumption that the glacier maintains a similar shape as it responds to mass-balance variations with changes in its area and volume. Timescales for redistribution of ice volume to maintain the characteristic shape of a glacier are expected to be much shorter than the response time of the glacier to massbalance changes (Jóhannesson et al., 1989;Harrison et al., 2001). This approach has its limitations for surge-type glaciers but may be expected to provide reasonable volume-change estimates over long time periods with substantial changes in volume, even for ice caps with many outlet glaciers that may surge at irregular intervals. In case the glacier surface and subglacial topographies are well known at one point in time, the computation of volume changes can be based on an accurate estimate of the glacier volume at that time and the method mainly relies on the assumption that there is a statistical relationship between volume and area changes (see Figure 4). Thus, uncertainty about the overall magnitude of the volume of the glacier in question does not affect the accuracy of the estimated volume changes in this case. We note that over the time periods considered in this article, repeated surface mapping and surface reconstructions of the glaciers, where available, have shown elevation changes that are small in the interior of the glacier and amplified toward the ice margin, as expected if the glacier maintains a geometrically similar shape as it adjusts toward a new geometry during variations in the climate (Hannesdóttir et al., 2015b;Thorsteinsson et al., 2017). FIGURE 4 | Volume-area scatter plot for the three largest ice caps in Iceland at various times since 1890. The blue markers are for Langjökull, green markers for Hofsjökull, and the black for Vatnajökull. The dotted gray line shows the least square fit of Eq. 1 through all volume-area points for these ice caps. Inset figures show a zoom in of the plot for the clusters with data from Hofsjökull and Langjökull (upper left corner) and Vatnajökull (lower right corner). Labels by the filled circles indicate dates when both area (Hannesdóttir et al., 2020) and volume are well established, while labels beside stars indicate dates of well-known area, with volume estimated from the deduced volume-area relation (dotted gray line). The 1890 values are obtained using area based on geomorphological evidence of the Little Ice Age maximum extent (Hannesdóttir et al., 2020). The 1890* and 1890** for Vatnajökull show the volume estimates with the reduced area to compensate for the large portion of surging outlets from Vatnajökull (see main text).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 The volumes of Hofsjökull and Vatnajökull in 1945 and 1890 and for Langjökull in 1890 were calculated based on the derived fit (shown with stars in Figure 4). The area of the glaciers in ∼1890 is based on geomorphological evidence of the maximum LIA extent in Iceland (Hannesdóttir et al., 2020). About half of Vatnajökull's glacier margin is substantially affected by surges (Björnsson et al., 2003). The geomorphological evidence indicating the maximum LIA extent in front of these surge-type outlets is due to an advance during a surge that increased the area of the glacier without much effect on the glacier volume. Using this area as input for the volume-area relationship (Eq. 1) is likely to result in an overestimate of the glacier volume in ∼1890 for two reasons. First, Vatnajökull did not at any one time reach the determined maximum LIA extent because the surges responsible for the geomorphological evidence of maximum extent did not occur simultaneously. Second, right after a surge, the glacier spreads over a larger area than before without any increase in volume, resulting in a thinner glacier with smaller volume than indicated by the volume-area relationship (Eq. 1). To take this into account, we reduced the area so it represents the likely ice-cap area when the surge-type outlets are on average near their mid-quiescent-period size. This corresponds to ∼500 m retreat from the determined LIA extent of Vatnajökull's surge-type outlets. Most of these advance between a few hundred meters and several kilometers during surges (Björnsson et al., 2003), the average probably being close to 1 km. The volume-area point marked 1890* for Vatnajökull in Figure 4 includes an area correction that corresponds to a 500 m retreat (area reduction by 100 km 2 ), and the point marked 1890** includes double this area correction (the point marked 1890 corresponds to data that have not been adjusted to reflect the impact of the surges on the area). The effect of surges on Hofsjökull and Langjökull is much smaller than on Vatnajökull and therefore not taken into account in this study.
From the estimated ∼1890 volume for the three ice caps and the 1945 volume for Vatnajökull and Hofsjökull, it is possible to estimate the mean specific mass balance for the periods between the volume estimates. The determined values for mean specific mass balance for two periods for Hofsjökull and Vatnajökull (1890/91 to 1944/45 and 1945/46 to 1969

Uncertainties of the Presented Mass-Balance Records
It is not straightforward to estimate the uncertainty of the massbalance estimates derived from the different observations and methods described above. A thorough consideration results in an  (Pálsson et al., 2012). A single value of 0.10 m w.e. a −1 is adopted for the smaller glaciers in 1945/46 to 2016/17. This corresponds to uncertainty values typically given for 10-20 year periods , while for periods exceeding 30 years, the uncertainty estimates are typically on the order of a few cm w.e. a −1 . The above uncertainties should be considered as possible biases and applicable for longer periods, exceeding decades. They neither reflect random annual errors in those records, nor annual deviation from long-term means for the geodetic and volume-area scaling results. When calculating the uncertainty of each time series shown in Figure 3G, the uncertainty for each contribution ΔC V , ΔC L , ΔC H , and ΔC O (corresponding to Vatnajökull, Langjökull, Hofsjökull, and "others", respectively) is derived by cumulating the assigned annual uncertainties. When calculating the uncertainty of the mass change of all Icelandic glaciers for the four IPCC periods (shown as horizontal lines in Figure 5), the uncertainties of the different contributions are considered independent. The uncertainty of the total mass change for 1970/71 to 2017/18, 1992/93 to 2017/18, and 2005/06 to 2017/18 is therefore For 1900/01 to 1989/90, the total uncertainty is due to how C O is calculated in 1890/91 to 1944/45 with the ratio F (as described above). Below we justify the above uncertainties and this approach.
In previous studies, we have assumed the uncertainty for a point measurement of the surface mass balance to be ∼0.30 m w.e, while the uncertainty of the specific annual mass-balance values is expected to be smaller due to the good coverage of the survey sites and the number of measurement stakes . The surface mass-balance record obtained with the glaciological method for Langjökull (in 1997/98-2003/04) has been compared with volume changes derived from geodetic mass-balance estimates (Pálsson et al., 2012). This comparison revealed ∼0.05 m w.e. a −1 higher mass loss from the glaciological method compared with the geodetic method. When adding the non-surface mass-balance component from Jóhannesson et al. (2020) of 0.055 m w.e. a −1 , the difference becomes ∼0.10 m w.e. a −1 . It is unfortunate that the more recent DEM [from 2004 used in Pálsson et al. (2012)] was obtained in mid-August so the surface melting until the end of the melt season (late September) was not accounted for. Therefore, the observed difference may partly be due to the lack of seasonal correction (that could be 5-20% of typical summer melt) when estimating the geodetic mass balance.
Thorough validation of the glaciological mass-balance record for Vatnajökull has not yet been carried out, but comparison between geodetic and glaciological surface mass balance for limited periods and areas has suggested differences that are comparable to the non-surface melt included here (Zamolo, 2019). Based on these comparsions, the uncertainties for the decadal averages for specific mass balance of Vatnajökull during the observed and modeled periods is determined to be 0.1 m w.e. a −1 , while the corresponding number for Hofsjökull and Langjökull is determined to be 0.07 m w.e. a −1 . These values are supposed to reflect the uncertainty of the average massbalance rate over a decade or longer period. For individual years, the uncertainty is assumed to be double this value, due to errors changing randomly from year to year.
The uncertainty of 0.1 m w.e. a −1 for the long periods obtained with volume-area scaling is supported by Figure 4. This value includes both the uncertainty of the volume used as input for the volume-area scaling and the uncertainty of the output values (shown with stars in Figure 4). This specific mass-balance uncertainty results in volume uncertainty of 127 km 3 and 72 km 3 for the Vatnajökull volumes in 1890 and 1945, respectively, when the uncertainties are cumulated from the date of the surface DEM (autumn 2010) used to calculate the volume (see Section 2.2). The corresponding minimum and maximum volume estimates for Vatnajökull in 1945 and 1890 (shown with error bars in Figure 4) are larger than the volumes estimated for the LIA maximum area (1890 star in Figure 4) and the doubled area correction due to surges (1890 ** star in Figure 4). We therefore consider the uncertainties of the specific mass balance for the period of volume-area scaling as a generous estimate. The same applies for the uncertainty of the specific mass balance for Hofsjökull and Langjökull.
The choice of the year 1890 as the time of maximum LIA extent of glaciers in Iceland for our analysis leads to some uncertainty in the mass-balance estimate for the first period after 1890, as the glaciers in fact reached the LIA maximum extent at different times. Most glaciers in Iceland reached their greatest historical extent during the LIA, with a maximum recorded in the late 19th century, although some glaciers reached a similar extent already during the 18th century (e.g., Thorarinsson, 1943;Geirsdóttir et al., 2009;Björnsson, 2017;Hannesdóttir et al., 2020, and references therein). Many glaciers started retreating from an advanced position near their LIA terminal moraines in the last decades of the 19th century, even if they reached the absolute maximum extent somewhat earlier. The volume of the glaciers will therefore have been close to the LIA maximum near the end of the 19th century, which we choose to represent with the year 1890, even for glaciers that reached their maximum earlier. We do not quantify the contribution of this uncertainty about the timing of the LIA maximum to the uncertainty of our mass-balance estimates as this can only be done on a glacier-by-glacier basis.
The mass change record of each glacier is constructed from three to four methods. Each method may have a constant bias, of a similar magnitude to the estimated uncertainties; the probability of the minimum mass change occurring for all periods (or alternatively maximum mass change for all Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 8 periods) is, however, smaller. By cumulating the uncertainties independent of the method, we assume this is possible. This generous uncertainty estimate is applied because of the partial dependency between methods, for example, the input data for the volume-area scaling is from both geodetic and glaciological data. Treating the methods as independent would therefore lead to underestimation of the uncertainties. The uncertainties for each of the contributions, ΔC V , ΔC L , ΔC H , and ΔC O , are assumed to be independent and therefore ΔC VLHO is calculated as the square root of a quadratic sum (Eqs 2 and 3). The more unrealistic assumption of full dependence between uncertainties, which would yield a direct sum of the uncertainties instead, would result in only ∼25% higher uncertainties. This is because ΔC V is by far the largest single contributor to the total uncertainty.

RESULTS
The mass change record of all glaciers in Iceland is shown in Figure 5. Before the glaciological year 1980/81, the observations do not allow the estimation of annual or decadal variability. The glaciological observations started on Hofsjökull in 1987/88 and annual variabilty in the period 1980/81 to 1987/88 is obtained from simulations from the HIRHAM5 snowpack model (Schmidt et al., 2019, see the Section 2). The annual variability is large: for the Vatnajökull record, the standard deviation of the observations is 0.75 m w.e. a −1 and for Langjökull and Hofsjökull 0.85 m w.e. a −1 and 0.78 m w.e. a −1 , respectively, during the period of glaciological observations on each glacier (see Figure 2).
Average mass change rates are computed for several selected periods (the reporting periods of the forthcoming IPCC AR6 assessment; see colored horizontal lines in Figure 5 Figure 3F shows that the cumulative mass change is primarily from Vatnajökull (365 ± 115 Gt) and that Langjökull and Hofsjökull have lost 63 ± 17 and 51 ± 13 Gt, respectively.
In the cold period 1980-1994 (see Figure 3E), 9 out of the 14 years show mass gain, but only 1 year after that (the massbalance year 2014/15). The high mass losses of 1996/97 and 2009/ 10 are conspicuous. The large mass loss in 1996/97 is due to the melting of ∼3.7 Gt of ice due to the subglacial Gjálp volcanic eruption in October 1996 (Guðmundsson et al., 2004), followed by a warm and sunny summer with low surface albedo due to dust precipitating onto the glacier surface. Both tephra and dust blown onto the surface of the glaciers enhance the surface melt due to lowering of the albedo (e.g., Möller et al., 2016;Wittmann et al., 2017). In the summer of 2010, a thin layer of volcanic tephra was distributed onto almost all Icelandic glaciers in the final phase of the Eyjafjallajökull eruption in April-May 2010 (Gudmundsson et al., 2012;Gascoin et al., 2017;Möller et al., 2019). Followed by an unusually warm and sunny summer, the tephra greatly enhanced melting Belart et al., 2019), especially in the accumulation zones of the glaciers, except on Eyjafjallajökull and large portions of Mýrdalsjökull, where the tephra layer became thick enough to insulate the snow and ice, reducing the melt rates (Dragosics et al., 2016).
The mass balance year 2014/15 was characterized by a long sequence of low-pressure systems arriving one after another through the winter, bringing large amounts of precipitation, followed by a cool summer with little melt, resulting in positive mass balance on all the glaciers. In the following three years, the mass loss rate was less negative than in the previous years, but the glaciological year 2018/19 was one of the most negative mass-balance years on record, due to the persistence of anticyclonic conditions during the summer of 2019 (Tedesco and Fettweis, 2020), which resulted in warm and sunny conditions from early spring. This led to early melting of a relatively thin winter snow layer, which in turn led to early exposure of lowalbedo regions in the ablation areas (Gunnarsson et al., 2020).

DISCUSSION
Although the glaciers in Iceland are located in a highly active volcanic region, they are useful monitors of climate variations. The detailed mass-balance record presented here is combined from glaciological observations, geodetic measurements, simulation with the HIRHAM5 snowpack model, estimates of non-surface mass balance, and results from an empirical volume-area scaling that are used to extend the record back to the time of maximum LIA extent of the glaciers as recorded by geomorphological evidence. The record spans 129 years, although the annual variability is not available until the last two decades of the 20th century. The record shows variability on decadal timescales with a period of near-zero mass balance in the 1980s and early 1990s before the onset of consistently negative mass balance on the order of −1 m a −1 that has prevailed since then. Close to half (−240 ± 20 Gt) of the total mass change (−540 ± 130 Gt) during the 129-year period occurred in 1994/ 95-2018/19, reflecting higher temperatures in this period (see Figure 3E) and is synchronous with glacier decline elsewhere in the world (Zemp et al., 2015). The record is valuable for studies of climate-glacier interactions and can be useful for validation of mass-balance models used to study glacier changes in other areas of the world where less data are available.
In the light of the limited geodetic observations ( Figures  3B,D) (Pálsson et al., 2012;Belart et al., 2020), the lengthchange observations back to 1930 (Eyþórsson, 1963;Sigurðsson et al., 2005), the glaciological measurements carried out on southeast outlets of Vatnajökull in the 1930s (Thorarinsson, 1940;Björnsson et al., 2013), modeling of the Vatnajökull SE outlet glacier Hoffellsjökull (Aðalgeirsdóttir et al., 2011), and the temperature evolution in Iceland during our study period ( Figure 3E), it is likely that a large proportion of the 20thcentury mass loss occurred during the ∼30 year period from the late 1920s to the late 1950s. For other time periods of the 20th Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 9 century, Icelandic glaciers were probably close to equilibrium on a decadal timescale. On those timescales, this study does not clearly indicate periods of substantially positive mass balance. In the period 1980/81-1993/94, which is probably the period with the most positive mass balance since the early 1920s, given the temperature history ( Figure 3E), the average mass gain only corresponds to 0.23±0.1 m w.e. a −1 for Vatnajökull, which is responsible for most of the 1.5 ± 1.0 Gt a −1 mass gain of Icelandic glaciers during that period.
Mass change records for Icelandic glaciers have been made from the glaciological observations provided to the World Glacier Monitoring Service (WGMS) database (Zemp et al., 2019) (a subset of the observations presented in Figure 3), from ICESat data (Nilsson et al., 2015)  is almost twice as large as our estimate, possibly due to signal leakage from mass changes of the neighboring Greenland Ice Sheet or an overcompensation for the effect of glacial isostatic rebound (e.g., Sørensen et al., 2017). A comparison of our results to the annual mass change rates of Zemp et al. (2019) and Wouters et al. (2019) is shown in Figure 6. This comparison shows that the interannual variability is generally well captured by both data sets, but some details are not; for example, the large ice melt due to the Gjálp eruption in October 1996 and the non-surface mass balance are not included by Zemp et al. (2019). The GRACE record (Wouters et al., 2019) has some years (e.g., 2006/07 and 2010/11) with more negative mass change, and others (e.g., 2005/06, 2011/ 12, and 2013/14) are less negative than our estimates, although the data points from our record are within the large uncertainty range of the GRACE values.
As discussed above, annual mass-balance data are not available prior to 1980/81, so we have used geodetic observations and the volume-area scaling to extend the record to the end of the 19th century. Zemp et al. (2019) extend the record for Icelandic glaciers back to 1950 using mass-balance observations from Storglaciären in Sweden and Storbreen in Norway. The obtained mass change rate for the period  in their record is about double the rate that we find here (−4.0 vs. −1.7 Gt a −1 ). Looking at decadal averages within this period, it appears that the difference is greatest in the period during which we assume near-zero balance, based on the temperature record ( Figure 3E), glacier length observations, and geodetic observations of the smaller glaciers . Our study thus shows that Scandinavian glaciers are not representative of glacier mass change in Iceland. Our results emphasize the importance of direct observations from glaciers located in Iceland for evaluation of global glacier variations.
The mass-balance record presented here includes the nonsurface mass balance due to geothermal melting, energy dissipation in the flow of water and ice, volcanic eruptions, and calving from the study of Jóhannesson et al. (2020). The previously estimated mass change rate of −9.5 ± 1.5 Gt a −1 for the period 1994period /95-2009period /10 (Björnsson et al., 2013, which included 0.5 Gt a −1 from geothermal melting (∼3% of typical ablation of the survey period), is less negative than the current estimate for the same period: −11.6 ± 0.8 Gt a −1 , now including the recent improved estimate of the contribution from other glaciers than the three largest, the non-surface melt, and calving in Jökulsárlón. As a fraction of the typical magnitude of the surface mass balance (∼−1 m w.e. a −1 on average for the glaciers in Iceland since 1995), the non-surface mass balance ranges from 5 to 38%, largest for FIGURE 6 | Comparison of mass change rates estimated in this study and studies based on glaciological observations provided to the WGMS database (Zemp et al., 2019;Zemp et al., 2020b) and the GRACE observations (Wouters et al., 2019), with the respective estimated uncertainties. The figure shows only the period when Vatnajökull, Hofsjökull, and Langjökull have all been monitored with glaciological observations. The contribution of other glaciers is from geodetic results , including an estimated annual variability (see Section 2).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 Mýrdalsjökull (31%) and the southern part of Vatnajökull (38%), and is therefore a non-negligible component of the mass balance.
Assuming that the geothermal melting and energy dissipation of water and ice flow were constant during the entire study period, the non-surface mass balance, excluding volcanic eruptions, contributed 1.5 Gt a −1 on average to glacier runoff during the period 1890/91 to 1994/95 and increased to 2.0 Gt a −1 on average after that (2.1 Gt a −1 if the Gjálp eruption in October 1996 is included). The time-averaged non-surface mass balance, therefore, explains ∼40% of the mass loss during the entire study period and ∼20% of the mass loss after 1994/95. There are strong similarities in the response of glaciers around the North Atlantic Ocean to atmospheric conditions, in particular at a multi-decadal time scale. This is illustrated by strongly negative mass balances and considerable retreat of glaciers in the 1930-1950s in Iceland (this study), on the east coast of Greenland (Bjørk et al., 2012), in Svalbard (Möller and Kohler, 2018) and for Hardangerjøkulen in Norway (Weber et al., 2019). More recently, the ocean around Iceland warmed after 1995 which correlates with the enhanced mass loss after 1995 in Iceland (Björnsson et al., 2013, this study) and Norway (Andreassen et al., 2020). The glaciers in Svalbard also shifted from a positive to a negative regime around 1990 (Østby et al., 2017;Schuler et al., 2020). The cooler oceanic conditions after 2010 (ICES, 2018) cooled the atmosphere and thereby reduced the mass loss in Iceland and Norway, and in Greenland the mass loss slowed down after a record mass-loss year in 2012 (Shepherd et al., 2019;Velicogna et al., 2020). The glaciological year 2018/19 shows one of the largest annual mass losses in our record. The Greenland Ice Sheet experienced a record mass loss in the same year (Sasgen et al., 2020) due to the persistence of anticyclonic conditions during the summer of 2019 (Tedesco and Fettweis, 2020).
With increased global temperatures, the summers in Iceland are warmer, resulting in longer ablation seasons, and in winter less precipitation falls as snow. With less snowfall on the glaciers, dirty ice appears earlier from beneath the snow in spring, which enhances the glacier ablation Gunnarsson et al., 2020). There is, however, a large variability in the melting due to several factors. In some years, the spring is cool, so glacier ice appears later from beneath the snow. There can be snowfall on glaciers during the summer months that rapidly raises albedo and reduces melt. Dust and tephra from volcanic eruptions are blown around in the highlands in dry periods and are often deposited on the glaciers, enhancing the melt (Wittmann et al., 2017).
With more detailed information about the past mass changes of Icelandic glaciers, models for projecting their future evolution can be improved. The mass loss from glaciers in Iceland has been projected to continue in the coming decades Marshall et al., 2005;Aðalgeirsdóttir et al., 2006;Guðmundsson et al., 2009;Aðalgeirsdóttir et al., 2011;Schmidt et al., 2019). Mass-balance models have been coupled with ice-flow models, and scenarios for future climate have been applied to simulate past and future mass changes. Hofsjökull and Langjökull, which currently have more negative specific mass balance than Vatnajökull and are both smaller in area and with less ice thickness, are likely to lose about 60 and 80% of their mass, respectively, until the year 2100 (Guðmundsson et al., 2009;Thorsteinsson et al., 2013). The larger ice cap Vatnajökull will survive longer; it is projected to lose 20-30% of its mass until the end of the century. Its future after that will depend on how much warming will be realized (Schmidt et al., 2019). These studies are based on mass-balance models that do not take non-surface mass balance into account and they therefore have a tendency to underestimate the future glacier decline. The projected mass losses toward the end of the 21st century are more rapid and persistent than the observations presented here. In the period 1890-2019, Vatnajökull, Langjökull, and Hofsjökull lost 12 ± 4, 29 ± 8, and 25 ± 6%, respectively, relative to their estimated mass in 1890. The future mass loss will both be due to the already realized temperature increase (Mernild et al., 2013;Vaughan et al., 2013;Marzeion et al., 2018) and the projected continued warming. The glaciers are now in disequilibrium with the already realized warming (Christian et al., 2018), because of their severaldecades-long response times (Jóhannesson et al., 1989;Harrison et al., 2001), and would continue to lose mass even if temperatures could be stabilized in the near future. The continued mass loss is dependent on future greenhouse gas emissions and whether, and if so how rapidly, they will be reduced.

CONCLUSION
Glaciers in Iceland are useful indicators of climate conditions in the middle of the North Atlantic Ocean. The results presented here add valuable information to global estimates of the response of glaciers to climate change in the past several decades. The study shows a total mass change of −540 ± 130 Gt (−4.2 ± 1.0 Gt a −1 on average) since the end of the LIA (∼1890), which corresponds to a 16 ± 4% loss of the LIA maximum ice mass. Almost half of the total mass loss has occurred since 1994/95 (240 ± 20 Gt, corresponding to 9.6 ± 0.8 Gt a −1 on average). The most rapid loss is observed in the period 1994/95 to 2009/10 (mass change rate −11.6 ± 0.8 Gt a −1 ). Since 2010, the mass loss rate has on average been ∼50% lower, with the exception of 2018/19, when one of the highest annual mass losses was observed (mass change rate −15.0 ± 1.6 Gt a −1 ). In 1930/31 to 1949/50, the average loss rates were probably close to the ones observed since 1994, while in 1980/81-1993/94, Icelandic glaciers had a period of small but significant surplus (1.5 ± 1.0 Gt a −1 ). For other periods of the study, the glaciers were either close to equilibrium or experiencing mild loss rates. This pattern of glacier evolution is similar across the North Atlantic region as shown by records from glaciers in Greenland (Bjørk et al., 2012), Svalbard (Möller and Kohler, 2018), and Norway (Weber et al., 2019;Andreassen et al., 2020).
After 1995, we have detailed glaciological observations, made on ∼60 survey sites on Vatnajökull and ∼25 survey sites on each of Hofsjökull and Langjökull, which show mass loss every year until 2014. In 2014/15, high winter precipitation and reduced melt during a short and cold summer caused a single anomalous year with positive mass balance. There is large interannual variability that often is impacted by volcanic eruptions enhancing the melt and dust or volcanic tephra blown onto the glacier surface from Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 523646 their sediment-rich vicinity. The glaciological year 2018/19 was among the most negative mass-balance years that were not significantly impacted by volcanic eruptions but by dust blown onto the glacier surface. The non-surface mass balance of glaciers in Iceland, as estimated by Jóhannesson et al. (2020), is non-negligible, accounting for about one-fifth of the mass loss since 1994. A part of this non-surface mass balance is caused by calving activity, which was insignificant in the first half of the 20th century, but has been gradually increasing with the ongoing retreat of the outlet glaciers located in over-deepened troughs . The calving will continue to increase as the glaciers retreat, and should, along with other non-surface mass-balance components, be taken into account in future projections of mass loss of glaciers in Iceland.

DATA AVAILABILITY STATEMENT
The regional mass-balance record constructed here is provided in Supplementary Tables S1, S2, together with mass-change rates for various periods. A subset of the mass-balance observations from Vatnajökull, Hofsjökull, and Langjökull is submitted to the WGMS database on annual basis. The geodetic mass balance estimates from Belart et al. (2020) are in the last version of the WGMS database (10.5904/wgms-fog-2020-08). The outlines of glaciers in Iceland published by Hannesdóttir et al. (2020) are submitted to the GLIMS database. Glacier length changes measured by members of the Icelandic Glaciological Society are available at spordakost.jorfi.is.

AUTHOR CONTRIBUTIONS
GA, FP, and EM designed the study, wrote the paper, and made the figures. HB and HHH designed, initiated, and led the Vatnajökull mass-balance work for decades and AG has facilitated the continuation of the mass-balance studies reported here. FP, EM, TT, TJ, AG, BE, HB, HH, HHH, and OS have conducted the in situ mass-balance measurements on the three main ice caps. HH, OS, EM, and JB contributed the area estimates at different times. JB, EM, EB, and TJ processed and contributed the geodetic estimates. LS contributed the simulations from HARMONIE-AROME. TJ contributed the non-surface mass-balance estimates. All authors contributed to discussions at various stages of the work and during the revisions of the manuscript.

FUNDING
This work and collection of data on which this work is based has been supported financially by and with the participation of the University of Iceland Research Fund; the National Power Company of Iceland; the Icelandic Public Road Administration; Iceland Glaciological Society; Reykjavík Energy Environmental and Energy Research Fund; three multinational European Union research projects TEMBA, ICEMASS, and SPICE; two projects supported by Nordic Energy Research: Climate and Energy (CE) and Climate and Energy Systems (CES); the Nordic Centre of Excellence SVALI (Stability and Variations of Arctic Land Ice), funded by the Nordic Top-level Research Initiative (TRI); and the Icelandic Ministry for the Environment and Natural Resources through the cooperative project Melting glaciers. SPOT6/7 data were obtained thanks to public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the program "Investissements d'Avenir" managed by the French National Research Agency. EB acknoweldges support from the French Space Agency (CNES).

ACKNOWLEDGMENTS
We express our gratitude to Þorsteinn Jónsson, Hlynur Skagfjörð Pálsson, Sveinbjörn Steinþórsson, Vilhjálmur Kjartansson, members of the Iceland Glaciological Society, and many others for their important contributions during decades of mass-balance measurements in the field. We thank three reviewers and the scientific editor Michael Zemp for constructive comments on the manuscript.