ARTICLES PUBLISHED ONLINE: 12 DECEMBER 2016 DOI: 10.1038/NGEO2863 Centennial glacier retreat as categorical evidence of regional climate change Gerard H. Roe1*, Marcia B. Baker1 and Florian Herla2 The near-global retreat of glaciers over the last century provides some of the most iconic imagery for communicating the reality of anthropogenic climate change to the public. Surprisingly, however, there has not been a quantitative foundation for attributing the retreats to climate change, except in the global aggregate. This gap, between public perception and scientific basis, is due to uncertainties in numerical modelling and the short length of glacier mass-balance records. Here we present a method for assessing individual glacier change based on the signal-to-noise ratio, a robust metric that is insensitive to uncertainties in glacier dynamics. Using only meteorological and glacier observations, and the characteristic decadal response time of glaciers, we demonstrate that observed retreats of individual glaciers represent some of the highest signal-to-noise ratios of climate change yet documented. Therefore, in many places, the centennial-scale retreat of the local glaciers does indeed constitute categorical evidence of climate change. A 2,000 1,000 0 Length (m) lpine glaciers are consequential and captivating elements of the Earth system that feature prominently in the lives of nearby communities1 . The nature of glacier motion was a research challenge for nineteenth-century physicists2,3 , and the late Holocene history of glacier margins has been a primary target of modern palaeoclimate science4 . Consequently, the centuryscale length history of several hundred glaciers is well known (for example, Fig. 1)4,5 . Whilst glacier mass balance (that is, areaaveraged accumulation minus ablation, ≡b (m yr−1 )) is a more direct measure of climate6,7 than glacier length, only a few dozen mass-balance records extend for more than two decades. The century-scale, near-global retreat of glacier fronts seems improbable without some coordinating global climate change. However, the formal statistical assessment of the role of climate change in glacier retreat has been limited to the numerical modelling of three individual glaciers, each with only a single set of model parameters8 ; and to a comparison of the global aggregate glacier mass loss in forced and unforced integrations of global climate models9 . By itself, any single glacier is a blunt statistical instrument. Each is a unique product of its local climate and landscape. Characteristic glacier-length response times of several decades10 imply only a few independent degrees of freedom in a centennial record, resulting in poor statistical resolving power to evaluate a trend. In part because of these factors, the most recent assessment of the Intergovernmental Panel on Climate Change (IPCC) concluded only it was ‘likely’ that a ‘substantial’ part of glacier retreat is due to anthropogenic climate change, a much weaker attribution than for other metrics of climate change11 . Here we introduce a method to combine glacier observations with the better-resolved local meteorological trends, which facilitates strong conclusions. The centennial-scale retreats of 37 widely dispersed glaciers have each necessarily required a climate change. And while the cause cannot be attributed purely from observations, the required climate changes are centennial-scale trends that are globally distributed. −1,000 −2,000 Hintereisferner, Austria Gries, Switzerland Nigardsbreen, Norway South Cascade, USA Koryto, Russia Franz Josef, New Zealand Frias, Argentina −3,000 −4,000 −5,000 1600 1650 1700 1750 1800 1850 1900 1950 2000 Year Figure 1 The global record of glacier lengths5 , for 158 glaciers with 20 or more individual observations (shown as dots). Coloured lines show the specific glaciers analysed in Fig. 4. The signal-to-noise ratio as a metric of glacier change We relate sL , the signal-to-noise ratio of glacier length (≡L), to sb , the signal-to-noise ratio of mass balance (≡b). Let 1L be the change in length over some period of time, and let σL be the standard deviation in the absence of any climate trend. Then sL ≡ 1L/σL . Likewise, sb ≡ 1b/σb . We first establish that sL is straightforwardly related to sb , and that the relationship is insensitive to uncertainties in glacier parameters. The result is robust and depends only on the fundamental property that glaciers integrate mass balance on timescales of a few decades. In a refinement of earlier models12,13 , previous work has shown that glacier flow on a sloping bed can be accurately emulated by a linear, third-order differential equation (Methods)14 . Let 1b(to ) be the change in mass balance due to a linear trend applied over a 1 Department of Earth and Space Sciences, University of Washington, Seattle, Washington 98195, USA. 2 Institute of Atmospheric and Cryospheric Sciences, University of Innsbruck, A-6020 Innsbruck, Austria. *e-mail: gerard@ess.washington.edu NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ 1 NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 ARTICLES Length response to trend 2 6 0 −500 −1,000 τ = 15 yr 4 150 3 5 4 1900 1950 50 2000 2 20 d 2 0 −2 1600 1700 1800 40 80 100 1900 2 0 −2 2000 1600 1700 Year e 60 Glacier response time (yr) Temp. (°C) Accum. (m yr−1) Year c 1 3 τ = 60 yr −3,000 1850 2 100 τ = 30 yr −2,500 5 6 4 Length (m) 7 5 0 −2,000 Amplification factor, γ 200 1 T (°C) −1,500 b Duration of trend (yr) a 1800 1900 2000 Year Length (m) 1,000 0 −1,000 No trend −2,000 1550 1600 Trend starting in 1880 1650 1700 1750 1800 1850 1900 1950 2000 Year Figure 2 The response of idealized glaciers to climate. a, The response to a warming trend (dashed line, right axis) of three idealized glaciers with response times of τ = 15, 30 and 60 yr (coloured lines, left axis). b, The amplification factor γ (τ , to ), in the relationship sL = γ · sb . The dashed red box covers the approximate range applicable to alpine glaciers. Thus, γ ' 4 to 6 for a wide range of the relevant parameter space. c, A 500-yr segment of synthetic, random white-noise accumulation (σP = 0.75 m yr−1 ), shading denotes 1, 2 and 3σ ranges. d, As for c, but for melt-season temperature (σT = 0.75 ◦ C), and with an imposed 1 ◦ C-per-century warming trend beginning in 1880. The blue line and shading have no warming trend. e, Response of a glacier with τ = 30 yr, β = 180. Due to the applied warming trend from 1880 to 2010, 1L = −1,700 m, σL = 460 m, giving sL = −3.7. In these simulations, anomalies in temperature and mass balance are related via b0 = −µT 0 , with µ = 0.65 m yr−1 ◦ C−1 . period to . Our model yields 1L(to ) = φ(to , τ ) · β · 1b(to ), where τ and β are functions of glacier geometry. The glacier response time τ = −H /bt , where H is a characteristic ice thickness and bt is the (negative) mass balance at the terminus12 ; β is the ratio of the glacier area to the product of H and the width across the terminus. The function φ(to , τ ) embodies the ice dynamics, and captures three distinct stages of adjustment: changes in interior ice thickness drive changes in ice flux at the terminus that, in turn, drive changes in glacier length12,14 . τ is central to a glacier’s response: Fig. 2a shows 1L(t) for a warming trend of 1 ◦ C per century, for three glaciers with different τ (and fixed β). Physically, τ controls how quickly a glacier responds to climate variations and also how strongly the glacier is restored to equilibrium12–14 : small-τ glaciers respond quickly but are less sensitive, whereas large-τ glaciers respond slowly but are ultimately more sensitive. These fundamental trade-offs are independent of the model used (see Supplementary Information), and mean that a century or so after a climate trend commences, the amount of retreat is relatively insensitive to several-fold variations in τ (Fig. 2a). Our model, and all equivalent models, also yield σL ∝ β · σb (Methods). Thus, sL can be written as: sL = γ (to , τ ) · sb (1) γ (to , τ ) is an amplification factor that depends only on the duration of the trend and the glacier response time. The counteracting 2 tendencies of initial responsiveness versus long-term sensitivity mean that γ (to , τ ) is quite insensitive to τ . For example, for a 130-yr trend, glacier length is a fourfold to sixfold amplifier of the massbalance signal-to-noise ratio, across a wide range of response times (Fig. 2b). The relative constancy of γ is key to estimating sL from meteorological observations. Note both 1L and σL are proportional to the parameter β and so it drops out of sL . The relationship between sL and sb can be illustrated with a synthetic example for typical climate trends and variability. Assume a 1 ◦ C-per-century increase in melt-season (June–September) temperature (≡T ) beginning in 1880, no trend in annualmean precipitation (≡P), and white-noise interannual variability (consistent with observations7 ) characterized by σT = 0.75 ◦ C and σP = 0.75 m yr−1 , respectively (Fig. 2c,d). A simple massbalance model is b0 = P 0 − µT 0 , where µ is a melt factor (= 0.65 m yr−1 K−1 )15 , and primes denote fluctuations about the long-term mean. After 130 yr of the imposed trend, sb = −0.65. For τ = 30 yr and β = 180, we get 1L = −1,700 m and σL = 460 m (Fig. 2e). Thus, the retreat is approximately three-and-a-half standard deviations (sL = 1L/σL ' −3.5), consistent with an amplification factor of γ ' 5.7 (equation (1) and Fig. 2b). An application to Hintereisferner, Austria We next present the steps of our analysis for Hintereisferner in the Austrian Alps (Fig. 1). Applying least-squares regression for the NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 TJJAS (°C) 4 2 0 −2 1880 100 Precip. (cm) b c ΔT = 1.4 °C; σ T = 0.7 °C 1900 1920 1940 1960 1980 2000 50 ΔP = −0.09 m yr−1; σ P = 0.13 m yr−1 0 1880 1900 1920 1940 1960 sb = −2.1(−2.8, −1.5) 1.0 0.5 0.0 −4 1980 2000 Year d 1.5 Probability density, h(sb) a ARTICLES e 6 −2 0 Mass balance signal-to-noise, sb 0.3 hs from Lobs Probability density Amplification factor, γ L 4 2 0 0 20 40 60 Glacier response time, τ (yr) L 0.2 hs combined L 0.1 0.0 −20 80 hs from P, Tobs, γ −15 −10 −5 Signal-to-noise ratio for length (ΔL/σ L) 0 Figure 3 Analysis for Hintereisferner (Austrian Alps, 46.8◦ N, 10.8◦ E). a, Melt-season (June–September) temperature from the Berkeley Earth data set of gridded station observations17 . The best-fit trend, the change since 1880 and the standard deviation are also shown. b, As for a, but from the Legates and Willmot data set of gridded annual-mean precipitation18 , with the trend extrapolated to 1880. c, The PDF of the signal-to-noise ratio for mass balance, sb (equation (2)), with median and 95% bounds also given. d, The blue shading on the x axis shows the PDF of uncertainty in glacier response time, τ ; the red curve shows the relationship between τ and the amplification factor, γ , for a 130-yr trend; the blue shading on the y axis shows the PDF for γ that results from the PDF for τ being projected onto the y axis via the red line. e, PDFs of length signal-to-noise ratio, sL , from different methods: blue, sL from observations of L; red, sL from T, P observations and the amplification factor γ (that is, equation (1)); orange, combined PDF using Bayes’ theorem. period 1880–2010, the glacier retreated 2,800 m with a standard deviation of 130 m about that trend. However, there are only three effective degrees of freedom in the record (Methods), so neither 1L nor σL is well constrained. Consequently, the probability density function (PDF) of their ratio (≡hsL Lobs ) is very broad (Fig. 3e and Methods): while it is extremely unlikely (<1%) that sL > 0, one cannot rule out the possibility its magnitude is very large (that is, sL < −20). In other words, there is not much information about sL from the glacier length alone. An independent approach to calculating sL uses equation (1). To that end, we build a simple mass-balance model using longduration, gridded instrumental observations16,17 of T and P, scaled by the observed variability of the much shorter winter (bw ) and summer (bs ) mass-balance records7,18 : b0 (t) = b0w (t) + b0s (t) = σbw T 0 (t) P 0 (t) + σ bs σP σT (2) where σ(x) is the standard deviation of x. Thus, the modelled mass balance matches the observed variance, and combines the observed signal-to-noise ratios of P and T . Normalizing by σP accounts for orographic-precipitation effects since fractional variations in precipitation are relatively uniform in mountainous areas. Between 1880 and 2010, Hintereisferner experienced strong warming: 1T = +1.4 ◦ C, with σT = 0.7 ◦ C (Fig. 3a); and some drying: 1P = −0.09 m yr−1 , with σP = 0.13 m yr−1 (Fig. 3b, extrapolated from a 100-yr record). Observations from the adjacent Vernagtferner18 give σbw = 0.22 m yr−1 , σbs = 0.42 m yr−1 . Thus, from equation (2), the median (and 95% bounds) for sb = −2.1(−2.8, −1.5) (Fig. 3c and Methods): sb is negative but with some uncertainty. Although the detrended b0 (t) is consistent with white noise, we evaluate the potential impact of climatic persistence in the Supplementary Information. For Hintereisferner19 , H ' 170 m and bt ' −6 m yr−1 , giving τ ' 30 yr. Formulae for τ vary in the literature12,13 , so we allow broad uncertainty; we represent its PDF by a gamma distribution with στ = τ/4. Figure 3d shows this nonetheless results in a very narrow PDF of γ , centred on γ ' 5.6. Hintereisferner thus acts as a near-optimal amplifier of the climate signal-to-noise ratio. From equation (1) the PDFs of γ and sb can be combined to give a second PDF for sL (≡hsL T ,Pobs ), based only on instrumental observations and the approximately several-decade response time of this glacier (Fig. 3e). This PDF rules out extremely negative sL (for example, sL ≤ −20) as inconsistent with observed climate trends (that is, sb ' −2) and the roughly sixfold amplification by the glacier length. These two PDFs of sL , one based on observations of L, and one based on γ and sb , are independent. Therefore, they can be combined using Bayes’ theorem (Methods)20 . The resulting median (and 95% bounds) is sL = −13(−17, −10) (Fig. 3e). This is an extraordinarily large magnitude compared with other documented climate metrics. For comparison, the local 1T near Hintereisferner is 2σ , and the global-mean, annual-mean 1T trend over the same period is 6σ . The result that 1L ≈ −13σL must not itself be directly equated to statistical significance because length variations are correlated in time. To proceed, we solve for the PDF of σL from the relation σL = 1L obs /sL , using the known observed retreat and the combined PDF for sL . For Hintereisferner, we find σL (and 95% bounds) = 220(170, 280) m. Although this estimate depends primarily on observations, it is consistent with calculations from modelling NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ 3 NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 ARTICLES h(ΔLnull) Hintereisferner, Austria σ L = 220 (170, 290) m sL = −13 (−10, −16) pLnull < 1% because of the observed drying in addition to the observed warming (Fig. 3a,b), and because the Bayesian step combining PDFs with length observations indicates a slightly more negative sL (Fig. 3e). ΔLobs An evaluation of glaciers worldwide Gries, Switzerland σ L = 260 (180, 420) m sL = −8 (−5, −12) pLnull < 1% Probability density Nigardsbreen, Norway σ L = 740 (470, 1,470) m sL = −5 (−2, −8) pLnull = 1% South Cascade, USA σ L = 610 (320, 2,740) m sL = −4 (−1, −7) pLnull = 6% Koryto, Russia σ L = 120 (80, 220) m sL = −7 (−4, −11) pLnull < 1% Franz Josef, New Zealand σ L = 460 (320, 720) m sL = −7 (−4, −10) pLnull < 1% Frias, Argentina σ L = 110 (80, 150) m sL = −15 (−12, −20) pLnull < 1% −4,500 −3,500 −2,500 −1,500 −500 500 Change in length (m) Figure 4 Analysis of glacier retreat from around the world (see Fig. 1). For each glacier, the PDF for 1L in any 130-yr period with no climate change is shown, and compared with the observed retreat represented by the vertical bar. Also given is the estimated standard deviation and signal-to-noise ratio (with 95% bounds), and the p value for the null hypothesis—the likelihood that the observed retreat occurred in the absence of a climate change. All PDFs have areas normalized to 1. See Supplementary Information for 30 other glaciers worldwide, and a full table of the analyses. (Methods). Finally, we test against a null hypothesis of no climate change. Using equation (1), the probability distribution of a given 1L occurring in any 130-yr period without climate change comes from combining the PDFs for each of the terms on the righthand side of: 1L null = σL · γ · sb null , where the PDF for sb null comes from the detrended mass-balance observations (Methods). We find the probability that the observed retreat comes from the null distribution is minuscule (pnull L = 0.001%, Fig. 4a); it is thus highly significant and must be attributed to a climate change. In general, the statistical significance of the glacier retreat may be larger or smaller than that of other local climate metrics. For example, although the global aggregate of the much shorter duration glacier massbalance observations is negative and significant at the 5% level7 , only 17/48 and 4/48 of the individual summer and winter mass-balance records longer than 10 years exhibit significant trends (all negative)7 . For Hintereisferner, sL and the statistical significance are so large 4 The foregoing analysis can be applied to any glacier with a long length record and known mass-balance variability. The results for seven widely distributed glaciers (Fig. 1) are shown in Fig. 4. Although the size of the retreats vary by more than a factor of five, for each glacier it is very unlikely that it could have occurred without climate forcing. The least significant retreat, pnull L = 6%, is South Cascade Glacier, which has a small median sb ' −0.7, resulting from a relatively small warming in a maritime climate where winter mass-balance variability exceeds that in summer7 (Supplementary Information). We have analysed a total of 37 glaciers worldwide, selecting those with the longest mass-balance and length records. All steps and results are detailed in the Methods and Supplementary Information. In several cases, nearby mass-balance records were used, and when length observations were too sparse to characterize the degrees of freedom, a uniform negative-definite hsL Lobs was assumed, consistent with the observed retreat (that is, 1L < 0) of all evaluated glaciers (Fig. 1 and Supplementary Information), and the fact that σL is positive. For the 37 glaciers, the median values of sL and σL range from −2 to −15 and from 120 m to 750 m, respectively. Our analyses represent an approximately tenfold increase in the number of glaciers for which σL has been estimated; three of our selected glaciers have previous model-based estimates of σL that lie within our observation-based ranges8,14 . Such estimates are valuable for putting reconstructions of past Holocene glacier fluctuations into context15 . We note that the response of glacier length to climate trends and to natural variability form complementary pairs: glaciers that are sensitive to climate change also tend to have higher variability. For 21 of the 37 glaciers, pnull ≤ 1%. Adopting IPCC L nomenclature, it is thus ‘virtually certain’ that the retreat of each of these individual glaciers required a climate change. A further seven have pnull L ≤ 5%. The least significant, Rabots Glacier in northern Sweden, has pnull L = 11%. Thus, for all but one glacier, it is ‘very likely’ (≥90%, IPCC) their retreat is attributable to climate change. These calculations and uncertainty estimates can undoubtedly be refined and improved by more sophisticated mass-balance models and by detailed numerical case studies with explicit valley geometry, and of course our analyses do not, on their own, speak to the cause of the required climate change. However, the decadal response time of glaciers means their centennial retreat is predominantly a response to the twentieth-century climate trends rather than being a dynamical recovery from any antecedent conditions, such as the putative Little Ice Age21 . The fundamental principle evinced here—that glaciers act as several-fold amplifiers of the signalto-noise ratio of local climate trends—is robust. The principle is not limited to glaciers. Any component of the climate system with a decadal timescale will damp high-frequency variations and integrate centennial-scale trends. However, glaciers are perhaps unique in combining a decadal timescale with a strong sensitivity and simple dynamical response to temperature, creating nearmaximum signal-to-noise ratios for centennial-scale climate change (Fig. 2b). In combination with climate and glacier observations, this has enabled our quantification of just how far individual glaciers have been driven from their pre-industrial states by climate change. Methods Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper. NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 Received 21 September 2016; accepted 18 November 2016; published online 12 December 2016 References 1. Nussbaumer, S. U. & Zumbühl, H. J. The little ice age history of the glacier des Bosons (Mont Blonc massif, France): a new high-resolution glacier length curve based on historical documents. Climatic Change 111, 301–334 (2012). 2. Forbes, J. D. Occasional Papers on the Theory of Glaciers (A & C Black, 1859). 3. Tyndall, J. & Huxley, T. H. Observations on Glaciers. Proc. R. Soc. Lond. 8, 331–338 (1857). 4. Solomina, O. N. et al. Holocene glacier fluctuations. Quat. Sci. Rev. 111, 9–34 (2015). 5. Leclercq, P. W. et al. A data set of world-wide glacier length fluctuations. Cryosphere 8, 659–672 (2014). 6. Braithwaite, R. J. & Zhang, Y. Relationships between interannual variability of glacier mass balance and climate. J. Glaciol. 45, 456–462 (1999). 7. Medwedeff, W. G. & Roe, G. H. Trends and variability in the global dataset of glacier mass balance. Clim. Dynam. http://dx.doi.org/10.1007/s00382-016-3253-x (2016). 8. Oerlemans, J. Holocene glacier fluctuations: is the current rate of retreat exceptional? Ann. Glaciol. 31, 39–44 (2000). 9. Marzeion, B., Cogley, J. G., Richter, K. & Parkes, D. Attribution of global glacier mass loss to anthropogenic and natural causes. Science 345, 919–921 (2014). 10. Leclercq, P. W. & Oerlemans, J. Global and hemispheric temperature reconstruction from glacier length fluctuations. Clim. Dynam. 38, 1065–1079 (2012). 11. Bindoff, N. L. et al. in Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) Ch. 10 (IPCC, Cambridge Univ. Press, 2013). 12. Jóhannesson, T., Raymond, C. F. & Waddington, E. D. Timescale for adjustments of glaciers to changes in mass balance. J. Glaciol. 35, 355–369 (1989). 13. Oerlemans, J. Glaciers and Climate Change (A. A. Balkema, 2000). 14. Roe, G. H. & Baker, M. B. Glacier response to climate perturbations: an accurate linear geometric model. J. Glaciol. 60, 670–684 (2014). ARTICLES 15. Anderson, L. S., Roe, G. H. & Anderson, R. S. The effects of interannual climate variability on paleoclimate estimates derived from glacial moraines. Geology 42, 55–58 (2014). 16. Rohde, R. et al. A new estimate of the average earth surface land temperature spanning 1753 to 2011. Geoinfor. Geostat. http://dx.doi.org/10.4172/gigs.1000101 (2013). 17. Legates, D. R. & Willmott, C. J. Mean seasonal and spatial variability in gauge corrected global precipitation. Int. J. Climatol. 10, 111–127 (1990). 18. Fluctuations of Glaciers Database (World Glacier Monitoring Service (2014); http://dx.doi.org/10.5904/wgms-fog-2014-09 19. Greuell, W. Hintereisferner, Austria: mass-balance reconstruction and numerical modelling of the historical length variations. J. Glaciol. 38, 233–244 (1992). 20. Annan, J. D. & Hargreaves, J. C. Using multiple observationally-based constraints to estimate climate sensitivity. Geophys. Res. Lett. 33, L06704 (2006). 21. Matthews, J. A. & Briffa, K. R. The little ice age: re-evaluation of an evolving concept. Geogr. Ann. 87 A, 17–36 (2005). Acknowledgements We are grateful to P. Green, K. Armour, D. Battisti and E. Steig for valuable comments and conversations. F.H. thanks the Institute of Atmospheric and Cryospheric Sciences, University of Innsbruck for financial support. Author contributions G.H.R., M.B.B. and F.H. planned the analyses, which G.H.R. performed. All authors contributed to the interpretation of the results and to writing the manuscript. Additional information Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to G.H.R. Competing financial interests The authors declare no competing financial interests. NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ 5 NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 ARTICLES Methods A three-stage model. A refinement of earlier analytical glacier models12,13 was developed in ref. 14, and showed that numerical models of ice flow on a constant slope could be accurately emulated by a third-order linear differential equation: d 1 + dt τ 3 L0 = β 0 b (t) 3τ 2 estimate sL from least-squares linear regression. Let hsL (sL ) Lobs be the probability density function (PDF) of sL , and let 1L and σL be the central estimates of 1L and σL , respectively. From ref. 22, and the standard formulae for uncertainties in regression parameters23 applied to glacier trends24 , it can be shown that: r (3) where L0 (t) is the length perturbation from some long-term, mean-state position, driven by mass-balance fluctuations, b0 (t). τ = −H /bt , where H is a characteristic thickness and bt is the (negative) mass balance at the terminus; β = Atot /(wH ), where Atot is the mean-state glacier area and w is the characteristic width at the √ terminus; and = 1/ 3. Reference 14 showed that equation (3) captures three essential stages of glacier adjustment: changes in interior ice thickness drive changes in ice flux at the terminus that, in turn, drive changes in glacier length. Equation (3) has analytic solutions that facilitate efficient analyses over a wide range of parameter space. hsL (sL ) Lobs = ν −1 F 12 (4) where φ(to , τ ) = τ r 3 −r (e − 1) + 1 + e−r +2 r 2 (5) and r = to / τ . We note that φ(to , τ ) → tτo (to − 3 τ ) in the limit that t0 τ . Equation (4) is used to calculate the curves in Fig. 2a. Glacier-length variance. Consider stochastic white-noise interannual fluctuations in mass balance (that is, random, year-to-year variability due to the vagaries of the weather, with no year-to-year persistence). Let the standard deviation of such fluctuations be σb . Such mass-balance fluctuations drive variability in glacier length, characterized by a standard deviation σL . Ref. 14 shows that: σL = β · ψ(τ ) · σb ψ(τ ) = τ (1 − κ)(1 + 4κ + κ ) (1 + κ)5 2 4 ν= 1t + 2 ! (9) n1t R∞ r(τ )dτ 0 (10) where n is the number of years in the record and 1t = 1 yr. The number of degrees of freedom in the length observations ranges from <1 up to ∼6. Low degrees of freedom yield broad distributions for hsL Lobs (for example, Fig. 3e), so that the most important information from the length observations is only that sL is not positive. When this method identifies ν ≤ 1, or when the length records are sparse, hsL Lobs cannot be calculated in this way. We instead stipulate a uniform negative-definite hsL Lobs , consistent with the observed retreat (that is, 1L < 0) of all glaciers analysed and the fact that σL is positive definite (see also Supplementary Information). Our second method for calculating sL (described below) uses climate observations with many more degrees of freedom, which produces a much narrower PDF (for example, Fig. 3e). When the two PDFs are combined, the narrower PDF dominates (Fig. 3e), and so our results and conclusions are not sensitive to the estimate of ν. Signal-to-noise ratio from climate observations and glacier amplification factor. The second method for calculating sL is to use the relationship from equation (8): sL = γ (to , τ ) · sb (7) b0 (t) = b0w (t) + b0s (t) = σbw Estimated variance for Hintereisferner. Taking Hintereisferner in the Austrian Alps as an example, reasonable values for parameters are19 : Atot ' 10.5 km2 , w ' 400 m, H ' 170 m; bt ' −6.5 m yr−1 , giving β ' 150 and τ ' 30 yr. With σb ' 0.5 m yr−1 from the adjacent Vernagtferner7,18 , equation (6) yields σL ' 230 m. This value is close to the central estimate of σL made from observed retreat of Hintereisferner and local climate trends: σL ' 220 m (Fig. 4). Signal-to-noise ratio, sL . Therefore, for a climate trend superposed on natural interannual variability, the signal-to-noise ratio for glacier length, sL , can be formed from equation (4) divided by equation (6): 1L φ(to , τ ) 1b(to ) = · = γ (to , τ ) · sb σL ψ(τ ) σb ν −1 12 (11) Calculating sb and its persistence. As described in the main text, the mass-balance model is: 21 with κ = 1 − (1t/ τ ) and 1t = 1 yr. We note that ψ(τ ) → (3τ 1t/16 )1/2 in the limit of τ 1t. sL = r Degrees of freedom in a record. We estimate the number of degrees of freedom, ν, in the glacier record following a procedure recommended in ref. 23. We fit a third-order autoregressive process to the length data14,23,25 and calculate the associated autocorrelation function, r(τ ). The effective degrees of freedom are then given by: (6) where ν −1 1L sL , ν − 1, 12 σL where F (x, ν − 1, µ) is a non-central t-distribution of the variable x with ν − 1 degrees of freedom, and non-centrality parameter µ (refs 22,24). The formula in equation (9) was verified in Monte Carlo simulations of time series with stipulated trends added to random realizations of noise. Response to a trend. Consider a linear trend in mass balance commencing at t = 0. Let 1b(to ) be the resulting change in mass balance after a time to has elapsed. The exact solution of equation (3) for the resulting length change, 1L(to ), is: 1L(to ) = φ(to , τ ) · β · 1b(to ) r (8) Contours of γ (to , τ ) are plotted in Fig. 2b. Signal-to-noise ratio from length observations. We considered glacier records from a compilation of 471 glacier-length histories worldwide5 . The data set uses Stineman interpolation in between individual length observations to produce annual time series, although our results are not sensitive to the interpolation method. We selected 37 glaciers in total (in each of five regions—see Supplementary Information), and focus on the period 1880 to 2010. The main criteria in selecting glaciers were long mass-balance records and a broad geographic distribution. However, we also preferentially selected glaciers with near-continuous length observations, although for some regions such as Asia and South America, this criterion could not always be met. The signal-to-noise ratio is given by sL = 1L/σL , where 1L is the change in length and σL is the standard deviation of the variability about that trend. We P 0 (t) T 0 (t) + σbs σP σT (12) Values of P 0 (annual mean) come from the Legates and Willmot data set17 . Values of T 0 (June to September in the Northern Hemisphere extratropics, December to March in Southern Hemisphere extratropics; annual mean in the tropics) come from the Berkeley Earth data set16 . The data were taken from the nearest land grid point to each glacier’s location. The data sets were chosen for their high spatial resolution, but alternative data sets show similar trends. The PDF of sb (≡hsb ) is calculated using the same non-central t-distribution as equation (9): r hsb (sb ) T ,Pobs = ν −1 F 12 r ν −1 1b sL , ν − 1, 12 σb r ν −1 12 ! (13) Here, the degrees of freedom are also calculated by fitting an autoregressive model to b0 (t). Only 3 out of the 37 glaciers analysed have a mass-balance time series that indicates any persistence (documented in the Supplementary Spreadsheet), with most being well characterized by white noise, consistent with mass-balance observations7,26 . Thus, for almost all of the records, the degrees of freedom are equal to the number of years in the record. The majority of the mass-balance trends are dominated by trends in summer mass balance due to the melt-season temperature trends, which is also seen in the much shorter records of actual glacier mass balance7 . Calculating the PDF of γ from the PDF of τ . The central estimate of τ comes from the relationship τ = −H /bt . For each glacier, the source references for the characteristic thickness H and terminus ablation rate bt are given in the Supplementary Spreadsheet. For the PDF of uncertainty in τ (≡hτ (τ )) we take a gamma distribution (ensuring τ is positive definite), with a standard deviation of τ/4 (Fig. 3d). Thus, the 95% uncertainty bounds on τ (that is, ±2σ ) are equal to τ itself—a broad and thus conservative estimate of the uncertainty. NATURE GEOSCIENCE www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ NATURE GEOSCIENCE DOI: 10.1038/NGEO2863 ARTICLES using hsL post (that is, equation (17)) and the fact that 1L is well constrained from observations. The PDF of γ (to , τ ) (≡hγ ) comes from equation (8), and applying the relationship: hγ (γ ) = hτ (τ ) ∂γ /∂τ (14) For our analyses to = 130 yr (1880 to 2010). The nature of γ is such that even very broad PDFs in τ map onto narrow PDFs of γ (Fig. 3d). The PDF of sL . Finally, the PDF of sL (≡hsL (sL ) T ,Pobs ) comes from combining hγ and hsb . The probability that sL lies between any given value sL and sL + dsL is given by: hsL (sL )dsL = X hγ (γ )dγ · hsb (sb )dsb (15) where the sum is taken over all possible pairs (sb , γ ) for which sb · γ = sL . Using dsb /dsL = 1/γ , and taking the limits that dsL , dγ , dsb → 0, the sum becomes an integral: hsL (sL ) T ,Pobs = ∞ Z hγ (γ ) · hsb (sL /γ ) · 0 1 · dγ γ (16) In other analyses in this study, the PDFs for all variables that are the product or ratio of two other variables are generated in a similar fashion. Where noted, the PDFs generated in this way have been verified with Monte Carlo simulations. Combining the PDFs for sL using Bayes’ theorem. hsL Lobs and hsL T ,Pobs are independent and so can be combined using Bayes’ theorem to yield an updated, posterior PDF for sL . We regard hsL Lobs as the prior PDF that is updated by new information from the observed sb . hsL T ,Pobs is equivalent to the likelihood function L(sL sb ), the likelihood of sL given the new information sb (ref. 20). For glaciers whose length history is too sparse (meaning if there are temporal gaps in the observations that are comparable to τ ) or if the degrees of freedom are less than one, we stipulate an uninformative, flat prior for hsL Lobs that is negative definite: hsL = constant for sL < 0, and hsL = 0 otherwise. This is consistent with the observed negative 1L of all glaciers analysed and the fact that σL is positive definite (see also Supplementary Information). Then the posterior PDF for sL is given by: hsL (sL ) post = k · hsL (sL ) Lobs · hsL (sL ) T ,Pobs (17) where k is a normalization constant independent of sL . An example of hsL (sL ) post for Hintereisferner is shown in Fig. 3e. Estimating σL from sL . We estimate the PDF of σL (≡hσL ) from the relationship: σL = 1L sL (18) Generating the PDF of the null hypothesis pL null . We evaluate the observed glacier retreat against a null hypothesis of no climate change. For no climate trend the PDF of sb (≡hsb null ) is governed by the standard Student’s t-distribution centred on zero24 . That is, equation (13) with 1b = 0. Thus, the PDF of the null hypothesis for sL , hsL null , is calculated from the relationship sL = γ · sb and combining hγ and hsL null . Finally, we have 1L = σL · sL , and so the PDF of the null hypothesis for 1L (≡h1L null ) comes from combining hσL and hsL null . We confirmed the validity of our calculations of h1L null by Monte Carlo methods: generating a PDF of 1Ls from 10,000 130-yr integrations of the three-stage glacier-length model driven by climate variability with no trend. Examples for h1L null are shown in Fig. 4 and in Supplementary Figs 3–8. Following the standard approach, we determine the statistical significance of a glacier’s retreat by calculating the probability, pL null , that a 1L could be consistent with the null hypothesis and yet lie beyond the observed retreat. Code availability. The three-stage glacier model code is available from: http://earthweb.ess.washington.edu/roe/GerardWeb/Home.html. All analyses were performed by implementing the equations presented in the paper, for the data sets documented above. Data availability. The following data sets were used in this study: glacier mass balance from: http://dx.doi.org/10.5904/wgms-fog-2014-09; glacier mass-balance analysis from: doi:10.1007/s00382-016-3253-x; glacier-length records from: http://www.the-cryosphere.net/8/659/2014/tc-8-659-2014-supplement.zip; Legates and Willmott precipitation data available from: http://www.esrl.noaa.gov/psd/data/gridded/data.UDel_AirT_Precip.html; Berkeley Earth temperature data available from: http://www.berkeleyearth.org/data. Other glacier parameters were taken from studies cited in the paper, and are documented in the Supplementary Information and in an accompanying spreadsheet. References 22. Nadarajah, S. & Kotz, S. Computation of signal-to-noise ratios. Commun. Math. Comput. Chem. 57, 105–110 (2007). 23. von Storch, H. & Zwiers, F. W. Statistical Analysis in Climate Research (Cambridge Univ. Press, 1999). 24. Roe, G. H. What do glaciers tell us about climate variability and climate change? J. Glaciol. 57, 567–578 (2011). 25. Schneider, T. & Neumaier, A. Algorithm 808: ARFIT—a Matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Trans. Math. Softw. 27, 58–65 (2001). 26. Burke, E. E. & Roe, G. H. The persistence of memory in the climatic forcing of glaciers. Clim. Dynam. 42, 1335–1346 (2014). NATURE GEOSCIENCE www.nature.com/naturegeoscience © ƐƎƏƖɥMacmillan Publishers LimitedƦɥ/ 13ɥ.$ɥ /1(-%#1ɥ 341#. All rights reservedƥ