Regional subdivisions<\/p>\n
We use the International Standardized Mountain Subdivision of the Alps (ISMSA)26<\/a> classification to define glacier regions in the European Alps. The ISMSA combines historical mountain range subdivisions of the adjacent countries and the division into Western and Eastern Alps, roughly divided by the Alpine Rhine Valley, Spl\u00fcgen Pass, and Lake Como. We calculate geodetic glacier mass-change rates for the two large divisions Western and Eastern Alps, and ten smaller glacierized subregions (see Fig. 1<\/a>), excluding regions with very small glacierized areas (<10\u2009km\u00b2).<\/p>\n The French Dauphin\u00e9 Alps (01 D-A) encompass the most western glacierized regions of the European Alps. The majority of glaciers are small (<0.5\u2009km\u00b2) and located close to the Massif des \u00c9crins. The Graian Alps (02 Gra-A) include glaciers in France, Italy, and Switzerland (>340\u2009km\u00b2). The largest glaciers are found at the Mont Blanc massif, surrounding the highest peak of the Alps (Mont Blanc, 4808\u2009m). The Pennine Alps (03 P-A) comprise the southern part of the canton of Valais (Switzerland) and the Aosta Valley (Italy). The second highest peak of the Alps (Dufourspitze, 4634\u2009m) is located in this region and also the second largest glacier coverage (>440\u2009km\u00b2). The Bernese Alps (04 B-A) are part of the canton of Bern and canton of Valais, Switzerland. The largest glacier of the Alps (Grosser Aletsch, 84\u2009km2 in 2000) is located within the Jungfrau-Aletsch mountain range with several summits above 4000\u2009m. The Bernese Alps also comprise the largest glacierized area in this study (>480\u2009km\u00b2). The Glarus Alps (05 Gla-A) and Lepontine Alps (06 L-A) border the Bernese Alps and Pennine Alps to the East, respectively. The glacierized areas are significantly smaller than in the adjacent regions (<60\u2009km\u00b2). The border between eastern Glarus and Lepontine Alps and Rhaetian Alps marks also the transition from Western to Eastern Alps. The Rhaetian Alps stretch across parts of Switzerland, Austria, and Italy. The majority of glacierized areas are located in the Bernina Range between Switzerland and Northern Italy (07 WR-A), \u00d6tztal, Austria (08 OR-A), and Stelvio National Park, Northern Italy (09 SR-A). The Tauern Alps (10 WT-A) are the most eastern glacier region and include the largest glacier of the Eastern Alps (Pasterze, 18\u2009km\u00b2 in 2000) and several other medium-sized glaciers surrounding the peaks of Gro\u00dfvenediger (3657\u2009m) and Gro\u00dfglockner (3798\u2009m). Average elevations and elevation ranges of glaciers in each subregion are summarized in Supplementary Table\u00a01<\/a>.<\/p>\n Area change calculation<\/p>\n Glacier outlines are calculated for three dates (2000, 2011, 2014) contemporanous with available DEM datasets. To avoid seasonal snow and extensive cloud coverage, late-summer (Aug\u2013Sep) Landsat images were selected, which represent the minimum glacier area at the end of the ablation period in the Alps. We use a total of 185 scenes from 1999\u20132001 (L5 TM & L7 ETM+), 2011 (L5 TM), and 2013\u20132015 (L8 OLI) to compute respective glacier areas corresponding to our DEMs. Band ratios (red\/shortwave-infrared)24<\/a> are created for each scene and converted to binary raster masks by applying manually selected thresholds. Thereafter, each raster is vectorized to create glacier polygons. To preserve comparability, ice divides between individual catchments are adopted from the Randolph Glacier Inventory V.6 (RGI)27<\/a> and RGI attributes such as glacier names and IDs are included.<\/p>\n For most timesteps and glaciers several summer acquisitions are available. To identify the most accurate glacier outline, all repeat coverage polygons are stacked individually for each glacier and polygons are selected automatically. Therefore, a set of criteria is evaluated for each glacier with several acquisitions available to derive the most reasonable glacier outline. As spatial parameters, the extent and area of the newly created outline is compared to the respective reference glacier statistics (of RGI) to identify exceptionally large, and thus unlikely, changes in glacier area and extent, e.g., due to partially cloud cover or seasonal snow patches. In addition, the Landsat Quality Assessment Band37<\/a> of the respective original image is used to derive a relative measure of so-called \u201cclear\u201d pixels within the direct vicinity of the glacier outline to further identify cloud or snow biased outlines. For glaciers with several available outlines of equal quality, the outline with the smallest respective glacier area is selected. In cases of incomplete spatial coverage of glacier catchments by all available acquisitions, polygon fragments from different acquisitions are combined and thereafter manually checked. In each case, the outline acquisition date(s) are preserved for each glacier. Outlines with a unique acquisition date are then used to calculate an area-weighted median acquisition date for each subregion, based on the glacier area enclosed by the respective outline, to derive regional area change rates during the observation periods.<\/p>\n The resulting glacier inventories are hereafter named according to their temporal composition. The first inventory (inventory 2000) includes primarily images from 2000 (61.3%) and 1999 (35.4%), the second inventory (inventory 2011) was created entirely from 2011 images and the third inventory (inventory 2014) mainly from 2014 (41.4%) and 2013 (35.0%) acquisitions.<\/p>\n Finally, each new inventory is visually inspected and remaining misclassified areas, such as patches of seasonal snow or proglacial lakes, are corrected manually based on false color composites, Landsat 7 panchromatic bands and elevation change raster. Debris-covered glacier outlines are also manually corrected according to respective elevation change fields and high-resolution satellite images (Google Earth). For the period 2013\u20132014, debris-covered glacier tongues are additionally compared with coherence estimates38<\/a> of Sentinel 1 image pairs from 2015 to distinguish debris-covered ice from rocks.<\/p>\n Elevation change calculation<\/p>\n We compute digital elevation models (DEM) and elevation change rates from Synthetic Aperture Radar data from the Shuttle Radar Topography Mission (SRTM) of the National Aeronautics and Space Administration (NASA) and of the TerraSAR-X add-on for Digital Elevation Measurement mission (TanDEM-X), operated by the German Aerospace Center (DLR) and Astrium Defense and Space. The SRTM dataset provides a consistent C-band DEM which was acquired during 11 days in February 2000 and covers all landmasses between 60\u00b0N and 56\u00b0S39<\/a>. We use the void-filled LP DAAC NASA Version 3 product with a ground resolution of 1\u2009arcsec40<\/a>. The TanDEM-X mission provides high-resolution X-Band acquisitions from 2010 onward with several complete coverages of the European Alps between 2011 and 2014. Within this study we produce two DEM mosaics from TanDEM-X CoSSC tiles for the periods 2011\u20132012 and 2013\u20132014 which are the only periods with enough acquisitions to cover the entire Alps. Whenever possible, we select TanDEM-X acquisitions from the same season as the SRTM DEM to minimize differences due to either radar signal penetration or snow accumulation at the acquisition time which can bias the elevation change measurement.<\/p>\n Each TanDEM-X elevation model is processed by using differential SAR interferometry, according to the workflow described in previous studies20<\/a>,21<\/a>,22<\/a>,23<\/a>. Differential interferograms are calculated using the void-filled SRTM DEM as reference surface. Subsequently, each interferogram is filtered and unwrapped by different algorithms (minimum cost flow & brach cut). The best results are selected manually and converted to elevation values by adding the reference surface heights. Thereafter, the newly created TanDEM-X DEMs are geocoded and coregistered to the reference DEM (SRTM DEM) to further reduce deviations between the DEM datasets. Each TanDEM-X DEM is vertically and horizontally coregistered to the reference DEM surface with an iterative process on stable areas. Those stable areas are selected by removing glacierized areas and slopes larger 15\u00b0. In addition, densely vegetated areas are excluded by using Landsat vegetation masks (Normalized Difference Vegetation Index). Subsequently, DEM mosaics are obtained from adjacent TanDEM-X DEMs and the acquisition date of each pixel is preserved alongside the elevation value. Finally, the coregistered TanDEM-X mosaic and the reference DEM (non-void-filled SRTM DEM) are differenced and change rates are computed using the respective TanDEM-X acquisition date and the date of the reference DEM. For the SRTM DEM we use the mean date (16-Feb-2000). Data voids due to gaps in the SRTM or TanDEM-X DEMs are filled by applying an elevation change versus altitude function, based on aggregated elevation change rates within 100\u2009m elevation bins41<\/a>. Contrasting to the hypsometric interpolation applied previously21<\/a>, we did not apply the three times the normalized median absolute deviation filter, which can introduce a bias on regional scales, particularly in the accumulation zones, as shown by a recent study42<\/a>. By testing of different filter approaches and manual inspection of the revealed results, a 1\u201399% quantile filter for each elevation bin was chosen instead to remove outliers in the region-wide hypsometric analysis. In addition, we remove steep slopes (>50\u00b0) where accumulation is negligible43<\/a>. As elevation reference, we use the void-filled SRTM DEM for the aggregation of elevation bins and hypsometric interpolation.<\/p>\n Geodetic mass change<\/p>\n Area and elevation change measurements are converted to mass budgets following the UNESCO definitions for glacier mass change estimates44<\/a>. We calculate geodetic mass changes for the periods 2000\u20132012 and 2000\u20132014, using the earliest glacier inventory (S1, 2000) as baseline and the inventories (S2) of 2011 and 2014, respectively, to determine a temporal mean glacier area (S3) for both observation periods as recommended by the UNESCO. Initially, elevation change rates are integrated for both periods over the respective maximum glacier area (Smax, spatial union of S1 and S2) and the change volume is calculated by multiplication of the derived elevation change rate and Smax. In addition, we apply a correction for SAR signal surface penetration (Vpen, see uncertainty section), which leads to an underestimation of volume change by the relative difference in signal penetration of the X- and C-band SAR. This bias volume due to signal penetration is then added to the measured volume change to derive the full glacier volume change. Subsequently, mass-change rates are estimated by applying a conversion factor assuming a mean density of 850\u2009\u00b1\u200960\u2009kg\u2009m\u22123 based on a study of alpine glaciers45<\/a>. We did not apply a variable density conversion as equilibrium line altitudes during the observation period are only available for a limited number of glaciers from a few regions in the Alps. Accumulation areas in many subregions are also very small and similar to the highest glacierized elevations. Therefore, we decided to use a constant conversion factor to provide a better comparability to other studies on geodetic mass change.<\/p>\n Finally, we determine specific-mass-change rates by dividing the volume changes by the respective temporal mean glacier area (S3) and multiplying with the density conversion factor.<\/p>\n Uncertainty assessment of geodetic mass change<\/p>\n We calculate the geodetic mass change uncertainty according to Eq. (1<\/a>) with \u0394M\/\u0394t being the mass change estimate, \u0394h\/\u0394t the average elevation change rate on the whole glacier surface, S1 and S2 the respective glacier areas at beginning and end of the observation period and p the applied volume to mass conversion factor:<\/p>\n $$\\delta _{\\Delta {\\it{M}}\/\\Delta {\\it{t}}} = \\sqrt {\\left( {\\frac{{\\Delta {\\it{M}}}}{{\\Delta {\\it{t}}}}} \\right)^2 \\times \\left\\{ {\\left( {\\frac{{\\delta _{\\Delta {\\it{h}}\/\\Delta t}}}{{\\frac{{\\Delta {\\it{h}}}}{{\\Delta {\\it{t}}}}}}} \\right)^2 + \\left( {\\frac{{\\delta _{{\\it{S}}_{_1}}}}{{{\\it{S}}_{_1}}}} \\right)^2 + \\left( {\\frac{{\\delta _{{\\it{S}}_2}}}{{{\\it{S}}_2}}} \\right)^2 + \\left( {\\frac{{\\delta _{\\it{p}}}}{{\\it{p}}}} \\right)^2} \\right\\} + \\left( {\\frac{{{\\it{V}}_{{\\mathrm{pen}}}}}{{\\Delta {\\it{t}}}} \\times {\\it{p}}} \\right)^2}$$<\/p>\n \n (1)\n <\/p>\n Which considers the following terms:<\/p>\n Error from the DEM differencing (including spatial autocorrelation and hypsometric gapfilling).<\/p>\n Error from glacier areas (independent errors from S1 and S2).<\/p>\n Uncertainty from volume to mass conversion using a fixed density.<\/p>\n Uncertainty from radar signal penetration.<\/p>\n The relative vertical precision of the elevation changes on the glacier surface and the hypsometric gapfilling contribute to the accuracy of the elevation change measurements (\u03b4\u0394h\/\u0394t). Elevation changes \u0394h\/\u0394t on stable areas, excluding glaciers, water and dense vegetation, are extracted and strongly deviating values are removed with a 2\u201398% quantile filter. Then, the elevation change values are aggregated in 5\u00b0 slope bins (Supplementary Fig.\u00a08<\/a>) and standard deviations (\u03c3\u0394h\/\u0394t) are computed to account for the dependance between surface slope and \u0394h\/\u0394t accuracy.<\/p>\n Each slope bin is filtered (2\u201398% quantile) to remove remaining artifacts. Eventually, the relative vertical precision of \u0394h\/\u0394t on the glacier areas is computed by weighting the obtained offsets for each slope bin by the slope distribution on glacier area (\u03c3\u0394h\/\u0394t AW) (Supplementary Table\u00a03<\/a>). To include the uncertainty contribution by spatial autocorrelation we generated semivariograms of 100,000 random \u0394h\/\u0394t samples on stable areas and derive a mean lag distance (dl) of ~312\u2009m. We follow a previous approach46<\/a> and estimate the accuracy of the region-wide average elevation changes according to Eq. (2<\/a>):<\/p>\n $${\\it{S}}_{{\\mathrm{cor}}} =\t \\, d_{\\it{l}}^2 \\times {\\it{\\uppi }}\\\\ \\delta _{\\Delta {\\it{h}}\/\\Delta {\\it{t}}} =\t \\, \\sqrt {\\frac{{{\\it{S}}_{{\\mathrm{cor}}}}}{{5 \\times {\\it{S}}_{\\mathrm{G}}}}} \\times \\sigma _{\\Delta {\\it{h}}\/\\Delta {\\it{t}}\\;{\\mathrm{AW}}} \\,\\,\\;{\\mathrm{for}}\\;{\\it{S}}_{\\mathrm{G}} \\, > \\, {\\it{S}}_{{\\mathrm{cor}}}\\\\ \\delta _{\\Delta {\\it{h}}\/\\Delta {\\it{t}}} =\t \\, \\sigma _{\\Delta {\\it{h}}\/\\Delta {\\it{t}}\\;{\\mathrm{AW}}}\\; \\qquad \\qquad \\,\\,\\,\\,\\,{\\mathrm{for}}\\;{\\it{S}}_{\\mathrm{G}} \\, < \\, {\\it{S}}_{{\\mathrm{cor}}}$$<\/p>\n \n (2)\n <\/p>\n where Scor is the correlation area and SG the mean glacier area (S3) multiplied by the empirical weighting factor 546<\/a>.<\/p>\n To account for errors due to misclassified glacier areas, we refer to a detailed comparison of automatically and manually classified outlines of alpine glaciers47<\/a>. The authors found a deviation in area of 3%, corresponding to a perimeter to area ratio of 5.03\u2009km\u22121. To account for different glacier geometries within our regional subdivisions (perimeter to area ratios ranging from 5.49 to 11.62\u2009km\u22121), we apply a scaling factor to better represent the larger area determination error in regions with very small glaciers. In addition, we estimate the accuracy of each glacier inventory individually to include both the uncertainties of the glacier area at the beginning (\u03b4S1) and at the end of the respective observation period (\u03b4S2) in our final error budget:<\/p>\n $$\\delta {\\it{S}}_1 = \\frac{{{\\it{r}}_{{\\it{P}}_1\/{\\it{S}}_1}}}{{{\\it{r}}_{{\\mathrm{P}}\/S_{{\\mathrm{Paul}}\\,{\\mathrm{et}}\\,{\\mathrm{al}}.}}}} \\times 0.03$$<\/p>\n \n (3a)\n <\/p>\n $$\\delta {\\it{S}}_2 = \\frac{{r_{{\\it{P}}_2\/{\\it{S}}_2}}}{{{\\it{r}}_{{\\mathrm{P}}\/{\\it{S}}_{{\\mathrm{Paul}}\\;{\\mathrm{et}}\\;{\\mathrm{al}}.}}}} \\times 0.03$$<\/p>\n \n (3b)\n <\/p>\n SAR signal penetration of the glacier surface strongly depends on the prevailing surface conditions during the SRTM and TanDEM-X DEM acquisitions. In general, surface penetration of X- and C-band microwave frequencies occurs more frequently during dry and frozen conditions whereas melting glacier ice leads to almost no penetration. When comparing X- and C-band, the signal penetration depth also differs due to different SAR frequencies. On cold dry ice, C-band frequencies can show maximum penetration depths of ~10\u2009m48<\/a>,49<\/a> while X-band penetration is smaller50<\/a>,51<\/a>,52<\/a>. Due to the unknown surface conditions and varying penetration depths of our DEM acquisitions we cannot quantify a specific penetration value for each DEM. Thus we estimate a bias volume (Vpen) due to the relative penetration difference of X- to C-band and correct the measured change volume over the periods 2000\u20132012 and 2000\u20132014 accordingly. We assume an increasing penetration bias from low to high altitudes. The SRTM DEM and most of the TanDEM-X scenes were acquired during winter months with temperatures well below 0\u2009\u00b0C at high altitudes and snow covered frozen glacier surfaces. At low altitudes, in the ablation zones, penetration depth was probably smaller because the surface ice was closer to the melting point53<\/a>. Therefore, we calculate a linear increase in signal penetration difference from 0\u2009m at minimum glacier elevation to 5\u2009m at maximum glacier elevation for each region as an approximation of the penetration depth difference between C- and X-band20<\/a>,21<\/a>. We use the 10% percentile of glaciated elevations, according to the glacierized area in 2000, as approximate minimum glacier elevation in each region. Eventually, we integrate Vpen in the error budget to account for a potentially under- or overestimated penetration bias in some regions.<\/p>\n Glacier-specific mass changes and comparison to glaciological measurements<\/p>\n In addition to the regional measurements, we calculate elevation changes for individual glaciers with glaciological mass balances for the periods 2000\u20132012 and 2000\u20132014. We aggregate elevation change rates within elevations bins of 10% of the glacier elevation range or 50\u2009m for glaciers with elevation ranges <500\u2009m or >500\u2009m, respectively41<\/a>. Outliers in each elevation bin are removed by applying a 1\u201399% quantile filter and data gaps in the hypsometric distribution of elevation changes are filled by a 3rd-order polynomial fit. To account for differences in signal penetration at the glacier surface, we apply on each glacier the regional correction function of the respective subregion. For the specific mass change, we use the temporal mean glacier area and a constant density of 850\u2009kg\u2009m\u22123 as for the regional mass changes. To estimate the uncertainty ranges of the individual glacier elevation changes, we use the respective regional elevation change uncertainties.<\/p>\n To compare our geodetic values with glaciological records, we calculate average annual mass changes of 25 glaciers29<\/a> with continuous measurements. For the glaciological averages we use values from the hydrological years 1999\/00\u20132010\/11 and 1999\/00\u20132012\/13.<\/p>\n Comparison of present mass-change rates to remaining ice volumes<\/p>\n