T he worldwide retreat of mountain glaciers is one of the most prominent impacts of climate change in public discourse and an established consequence of anthropogenic climate change [1][2][3] , with Andean glaciers among the fastest retreating 3,4 . Mountain glacier lengths have decadal-to-multidecadal response times 5,6 and act as lowpass filters of climate 6,7 . They therefore represent a purer signal of anthropogenic climate change than most other metrics. However, any individual glacier's response to climate change depends on its geographical and climatic setting 1,2,8 . An assessment of the human contribution to an individual glacier's retreat and hence changing glacial lake outburst flood (GLOF) hazard must quantify these factors and place the observed change in the context of the natural variability expected in the absence of anthropogenic climate change.

GLOFs are a dangerous hazard in deglaciating alpine environments, and downstream inundation from large outburst floods have caused serious damage to human settlements 9 and infrastructure 10 . GLOFs can occur from englacial, subglacial, supraglacial and proglacial lakes. GLOFs from proglacial lakes that form in the wake of a retreating glacier, such as Lake Palcacocha, are a response to the disequilibrium imposed on glacial landscapes by climate change 11,12 . They occur when an unstable moraine dam is breached or a moraine-overtopping wave is induced by a triggering event, such as a landslide or an avalanche 13 . However, the anthropogenic influence on a specific GLOF has hitherto not been evaluated 9 . In the mid twentieth century, new proglacial lakes formed throughout the Cordillera Blanca as glacier termini receded 14 and multiple fatal GLOFs were recorded, including those from Lakes Ayhuinyaraju, Jancururish, Tumarina and Palcacocha 15 . The ongoing GLOF hazard continues to require extensive mitigation and engineering efforts 16,17 . We focus here on Lake Palcacocha, which currently threatens the city of Huaraz with a GLOF 3,[11][12][13]16,18,19 .

Lake Palcacocha (9° 23′ 49″ S 77° 22′ 47″ W) sits at the foot of Palcaraju glacier. Lichenometric dating of Palcaraju glacier's terminal moraine indicates a seventeenth century formation, which sets the maximum extent of the glacier in the past few centuries 15,20 .

Photographic records of Palcaraju glacier (Fig. 1) show the 1939 terminus set well back from the terminal moraine, with an overspilling Lake Palcacocha between the glacier's terminus and terminal moraine. In 1941, the breaching of Lake Palcacocha's moraine dam resulted in the Cordillera Blanca's deadliest GLOF of the twentieth century 15 . The resulting debris flow destroyed one-third of Huaraz and killed at least 1,800 people 12,14,21 . A limited expansion of Lake Palcacocha was observed between 1948 and 1995, and the lake area increased from 0.06 to 0.08 km 2 . The rapid retreat of Palcaraju glacier since 1995 expanded the lake to an area of 0.49 km 2 in 2018 22 . Two concrete overflow spillways were installed in the terminal moraine between 1972 and 1974. Up to 2011, the lake surface elevation remained constant and the lake volume grew as Palcacocha expanded into the space previously occupied by the retreating glacier. In 2011, siphons were installed to lower the water surface by 3-5 m, although the GLOF hazard remains high 12 . The most recent bathymetric survey found that the volume of Lake Palcacocha was 17,403,353 m 3 in 2016 23 , a 34-fold increase from the volume of the post-GLOF lake in 1941, which makes it the fourth largest of the 35 glacial lakes with installed safety features in the Cordillera Blanca 16 .

Any evaluation of the contribution of anthropogenic greenhouse gas and aerosol emissions to a changing GLOF hazard must consider the whole causal chain. Our study has three parts. We evaluated: (1) the anthropogenic contribution to the observed temperature trend, (2) the role of this trend in the observed retreat of Palcaraju glacier (which comprises numerical modelling of the glacier and a statistical comparison of the observed retreat with the natural length variability) and (3) the role of the glacier's retreat in reshaping the valley and the resultant change in the GLOF hazard.

Globally, virtually all long-term temperature variations since the mid-nineteenth century have been attributed to external forcing 24 . Following the approach taken by Haustein et al. 25 A potential glacial lake outburst flood from Lake Palcacocha (Cordillera Blanca, Peru) threatens Huaraz, a city of 120,000 people. In 1941, an outburst flood destroyed one-third of the city and caused at least 1,800 fatalities. Since pre-industrial times, Lake Palcacocha has expanded due to the retreat of Palcaraju glacier. Here we used observations and numerical models to evaluate the anthropogenic contribution to the glacier's retreat and glacial lake outburst flood hazard. We found that the magnitude of human-induced warming equals between 85 and 105% (5-95% confidence interval) of the observed 1 °C warming since 1880 in this region. We conclude that it is virtually certain (>99% probability) that the retreat of Palcaraju glacier to the present day cannot be explained by natural variability alone, and that the retreat by 1941 represented an early impact of anthropogenic greenhouse gas emissions. Our central estimate is that the overall retreat is entirely attributable to the observed temperature trend, and that the resulting change in the geometry of the lake and valley has substantially increased the outburst flood hazard.

Way (HadCRUT4-CW) dataset 27,28 (Fig. 2a). The HadCRUT4-CW trend is well supported by regional observations (Methods). To quantify the anthropogenic contribution to regional temperature change, we regressed the observed temperatures in HadCRUT4-CW onto the Global Warming Index (GWI) 25 , an estimate of the anthropogenic contribution to global externally forced temperature change, over the period 1940-2018 (Methods). This yielded region-specific scaling factors for the GWI. Regional anthropogenically forced warming is the product of this scaling factor and the GWI anomaly relative to 1880. For 1989-2018, it was found to be equal to 95% (85-105%, 5-95% range, all percentages rounded to the nearest 5%) of the observed temperature change over the same period. Uncertainty bounds are given by the intermodel range of regressions of the relevant grid cell in the control runs (with non-evolving pre-industrial conditions) of the Coupled Model Intercomparison Project Phase 5 (CMIP5) against the GWI. Figure 2a shows the observed temperature change for the HadCRUT4-CW grid cell that contains Palcaraju glacier (black) and the temperature change in the same grid cell attributed to human greenhouse gas and aerosol emissions (orange). The correlation Fig. 2 | Regional observed and attributable temperature anomalies, and impacts on decadal mean mass-balance profile. a, Annual mean temperature (HadCRUT4-CW) expressed as anomalies from 1850 to 1880 and the resulting ELA, and the contribution of anthropogenic forcing to the observed temperature change, with its 5-95% uncertainty ranges (shaded orange). b, The resulting mass-balance profiles, plotted as decadal means, with the 1σ and 2σ mass-balance range around the base profile shown in grey. The glacier's 1880 profile (Fig. 3) is in mass balance for the base profile.

Nature GeoscieNce between regional temperatures and global anthropogenic warming with the estimated response to the natural external climate drivers (solar variability and volcanic activity) is negligible.

The climate that a glacier experiences can be characterized by the vertical profile of its mass balance (that is, accumulation minus ablation) 29 . We used a mass-balance profile (Fig. 2b) for Palcaraju glacier, which, for the current climate, matches that observed on Shallap glacier (9° 29′ 24″ S, 77° 20′ 24″ W, a well-observed glacier approximately 10 km away, with similar aspect to that of Palcaraju) 29 . The equilibrium line altitude (ELA, ~4,950 m in the current climate) denotes the altitude above which there is net accumulation and below which there is net ablation. A simple and robust method to quantify the impact of climate change on glacier mass balance is via geometrical shifts of the mass-balance profile in response to temperature and precipitation changes [30][31][32][33] . The climatological precipitation is important for the shape of the mass-balance profile but, on century timescales, precipitation trends are small and temperature changes dominate the mass-balance changes 32,34 (see also Methods and Supplementary Information). Temperature changes are implemented through a vertical shift of the ELA based on the temperature lapse rate (here, approximately 5.5 K km -1 , which gives an ELA sensitivity of approximately 180 m K -1 ) (refs. 32,35 ). Detailed surface-energy-balance models show a close agreement with this value for glaciers in the Cordillera Blanca 32,36 .

Figure 2b shows the decadal mean mass-balance profiles produced in response to the HadCRUT4-CW observations since 1880.

Tropical glaciers are characterized by steep vertical gradients in ablation, which renders them particularly sensitive to climatic perturbation 35 . Starting from an 1880 baseline, the point at which anthropogenic warming begins to emerge (Fig. 2a), warmings of 0.25 and 1.3 K (the respective anthropogenic contributions in 1941 and 2019) raise the ELA by 45 and 235 m, respectively. At an elevation of 4,600 m, characteristic of the early twentieth century ablation zone, these temperature changes increase the net ablation rate by 30 and 150%, respectively. Given the linear relationship between temperature and ELA, human-induced warming (as previously stated) and the change in ELA are equal to 95% of the observed temperature and modelled ELA changes, respectively.

Next, we used a numerical flowline model of ice dynamics (Methods) to evaluate the thickness profiles of the glacier during its retreat, the interaction of a rising ELA with the valley geometry and the natural glacier fluctuations in the absence of climate change. In the absence of direct observations, we estimated that the nineteenth century glacier terminus was situated 4.25 km along the flowline (Methods). Uncertainty in the 1880 terminus position has no effect on the present-day terminus position due to the geometry of the valley and the short response time of the modern glacier (Supplementary Information and Supplementary Fig. 1).

We modelled the response of Palcaraju glacier to regional annual temperature anomalies from three sources: HadCRUT4-CW (Fig. 2a), Berkeley Earth Land/Ocean observations 37 and the European Centre for Medium-Range Weather Forecasts' twentieth century reanalysis (ERA-20C) 38 .

For all three temperature time series, the numerical model gives an overall retreat similar to that observed (Fig. 3a,b), with the timing of the period of most rapid retreat being less predictable. The planform insets in Fig. 3a show the regions of the glacier with net accumulation in blue and net ablation in red, and highlight that this non-linear retreat is a geometrical consequence of the steep vertical gradient in mass balance and shallow slopes in the lower valley. Applied to the 1880 profile, a warming of 0.25 K causes the net ablation below the ELA to increase by 40%, from 1.8 m yr -1 to 2.5 m yr -1 , which cannot then be balanced by the 1.1 m yr -1 of net accumulation above the ELA. For 1.3 K of warming, the average net ablation rate increases by 230% to 5.9 m yr -1 , and there is a 10% decrease in net accumulation above the ELA to 1.0 m yr -1 . Rising temperatures have expelled the glacier from the lower valley.

In the observations and our simulations, the glacier retreats rapidly between lengths of approximately 3.0 and 2.7 km, taking four years in observations, and one, six and six years in our modelled response to the HadCRUT4-CW, Berkeley Earth and ERA-20C, respectively, time series. This period of rapid retreat is present in all the simulations, although its timing is sensitive to the warming history. It occurs in our model when thin ice above the slope break completely ablates to isolate a patch of stagnant ice on the lower slope. For the real glacier, the timing and nature of this retreat also depends on the glacier-lake interactions and other processes not represented in our simulations. Our numerical model lacks lakeice dynamics (which include glacier calving and thermal conduction) and debris cover (Methods), which probably affects the match between the transient evolution of the observed and modelled termini (Fig. 3b).

Could the observed retreat have happened without anthropogenic global climate change? Interannual climate variability will cause glacier fluctuations even in the absence of large-scale climate change 6,39 . The standard deviations (1σ) of regional annual mean temperature and precipitation are approximately 0.5 K and 0.2 m yr -1 , respectively, with some interannual persistence, probably due to the influence of El Niño (Methods). Applying this climatic variability to the 1880 glacier profile in the flowline model gives an interannual variability in ELA and mass balance of 100 m and 0.75 m yr -1 , respectively (both 1σ), consistent with observations of similar glaciers 40 . We generated 10 kyr of synthetic climate variability and found centennial-scale length fluctuations with a standard deviation (σ L ) of 230 m (Figs. 3b and4a). Figure 4b,c shows histograms of the trends over 60 and 140 years from the model output, compared with the observed retreat of Palcaraju between 1880 and 1940 or 2018, estimated at 550 and 1,850 m, respectively. In both cases the rates of observed retreat are larger than any of the model output, which indicates that a large-scale climate change is a necessary condition for the glacier's observed length change to occur.

We then quantified the change in glacier length that can occur from natural variability alone (ΔL| null ), which includes uncertainty in the glacier response time (τ) and the standard deviation of the glacier length (σ L ). This analysis accounts for glaciers' amplification of the signal-to-noise ratio of climate change 2 (Fig. 4d and equation (4); Methods) and generates probability density functions (PDFs) for trends over 60 and 140 years in the absence of climate change (Fig. 4e,f; Methods). Given our assumptions, we conclude that it is virtually certain (>99% probability) that the observed retreat of Palcaraju glacier could not have occurred due to natural variability alone and therefore that the observed large-scale climate warming that we attribute to human influence is a 'necessary cause' 41 of the observed retreat, both to 1940 and to the present.

The PDF of ΔL| null is symmetrical and centred over 0 m. Thus, over a 140 year period, the probability that the glacier would have extended is equal to the probability that it would have retreated, so our central (mean and median) estimate for ΔL| null is 0 m. Consequently, our central estimate for the attributable fraction of the observed retreat to the observed temperature trend is 100%.

Above, we show that anthropogenic warming equals 95% of the observed temperature change in 1989-2018.

Further, for the 1940 retreat not to have been the result of early anthropogenic warming (and be explained only by natural variability) at the 90% confidence level, the pre-anthropogenic equilibrium terminus must have lain within 200 m of the observed 1940 position (Fig. 4e and Supplementary Fig. 2d). This is inconsistent with the early twentieth century retreat observed elsewhere in the Cordillera Blanca and the glacier-length sensitivity to climate change required for the formation of the seventeenth century moraine (Methods, and refs. 14,[42][43][44] ). It would also require some unknown factor to counter the impact of the 0.25 K of anthropogenic warming on the mass balance (Fig. 3a inset).

Our results establish that the present geometries of the glacier and lake are the result of anthropogenic climate change. What are the implications for a GLOF hazard? A GLOF hazard is the product of the probability of a GLOF event occurring and the event magnitude 45 . Previous work has shown that, although a catastrophic breach of the moraine dam is unlikely, the present GLOF hazard from Lake Palcacocha is very high due to the substantial likelihood of an avalanche or landslide inducing a moraine-overtopping wave 12,19 . We used two ranking methodologies to evaluate changes in the GLOF hazard from Lake Palcacocha between the present day and the nineteenth century (Methods and Supplementary Tables 2 and3). The ranking schemes indicate that the greatest influences on the increasing GLOF hazard are the steepness of the glacier snout, the growing lake area and the glacier's retreat. The geometrical method of Wang et al. 46 was developed for settings in which a moraine-overtopping wave is the likely GLOF mechanism 46 . In line with previous applications of this methodology to Palcacocha 11 , we confirmed that the present GLOF hazard is categorized as 'very high' and found that it has increased from 'medium' in the nineteenth century. We compared these results with a geomorphological assessment methodology adapted from Bolch et al. 47 . This also ranks Lake Palcacocha in the highest GLOF hazard category, with the nineteenth century GLOF risk categorized as 'medium' . As a result of this increase in GLOF hazard, Lake Palcacocha is considered a serious threat to Huaraz, which compels the local authorities to implement hazard mitigation measures 19 .

The retreat of Palcaraju glacier, and the resultant expansion of Lake Palcacocha, has increased both the probability and potential magnitude of a GLOF. In addition to the criteria assessed by the two hazard assessment methodologies employed here, permafrost degradation in response to rising temperatures can decrease slope stability, which increases the probability of rockfalls and landslides into glacial lakes. Further, the retreat of glaciers, such as Palcaraju, to steep mountain headwalls brings proglacial lakes closer to icy mountain flanks that have decreased stability due to glacier de-buttressing 47,48 . The increase in lake size since the nineteenth century has also increased the likely GLOF volume, although neither methodology accounts for this change (the methodology of Bolch et al. 47 includes area, but the lake is already in the highest area class in the nineteenth century).

The GLOF hazard posed by Lake Palcacocha to Huaraz has increased substantially due to the lake's expansion as Palcaraju glacier retreats and has now become a critical threat to Huaraz, which requires engineered mitigation efforts. Given the assumptions in this three-step attribution, we found: (1) anthropogenic warming equals 95% (85-105%, 5-95% confidence interval) of the observed warming in 1989-2018, (2) our central estimate is that the retreat of Palcaraju glacier is entirely the result of this warming, (3) the retreat cannot be explained by natural climatic variability alone and (4) this retreat has substantially increased the GLOF hazard for Huaraz. Despite some uncertainty about the 1880 position of the glacier

Nature GeoscieNce terminus, our analysis also provides strong evidence of anthropogenic influence on the 1941 GLOF, a mid-twentieth century humanitarian catastrophe. For the attribution of the overall retreat of Palcaraju glacier, the nineteenth century position of the terminus is of lesser importance, and the present-day threat to Huaraz is a direct consequence of the anthropogenically driven retreat of Palcaraju glacier.

Our analyses demonstrate the utility of risk-based approaches to evaluate the anthropogenic contribution to slow-onset, persistent hazards. The set of analyses presented in this article could be applied to assess anthropogenic influence on the GLOF hazard in any deglaciating alpine environment. Application of these methods in other locations would provide regional and global insight into the human influence on present GLOF hazards and indicate how a GLOF hazard may evolve under future climate change and where adaptation measures may need to be implemented. 2 ). We accounted for uncertainty in the value of τ in 1880 in the PDF on the horizontal axis (grey shading), which gives rise to a PDF of values for γ for each of the time periods (shaded PDFs on the vertical axis). e, The PDF for ΔL| null in any 60 yr period in a climate with interannual temperature variability but no climate change trend, compared with the observed retreat between 1880 and 1940 (vertical bar). f, As e, but for a 140 yr period and as compared with 1880-2019.

Nature GeoscieNce

Temperature data and lack of precipitation trend. Given the minimal seasonal temperature variations and consequent year-round ablation from tropical glaciers 29,35 , we modelled the evolution of Palcaraju glacier in response to annual mean temperatures. We primarily used the UK MetOffice HadCRUT4-Cowtan/ Way (HadCRUT4-CW) dataset 27,28 as it provides continuous observations since 1850 and therefore a better pre-industrial baseline temperature for Peru than those of other observational datasets. Of the datasets considered, it also has an interannual temperature variability closest to the observed one. The standard deviation of annual mean temperature at the closest high-altitude meteorological station with an adequate record (the Mayor General FAP Armando Revoredo Iglesias Airport, 2,700 m) is 0.62 K compared with 0.64 K for HadCRUT4-CW over the same period (1963-2018).

The baseline temperature was taken as the mean of 1850-1880. The observed warming between the baseline and most recent 10 yr period (2009-2018) in HadCRUT4-CW is 1.01 K. Station observations in the Tropical Andes give temperature trends of +0.31 K per decade for 1969-1998 and +0.13 K per decade for 1983-2012 50 . Given these observed trends and the amplification of warming at higher altitudes 51 (the mean elevation of the HadCRUT4-CW grid cell is 1,333 m), it is likely that our estimate of regional warming in the HadCRUT4-CW dataset is conservative. These factors do not affect the proportion of the trend attributable to anthropogenic influence. Further, surface-temperature simulations from the Hadley Centre regional climate model version 3P (HadRM3P) coupled to the global atmosphere model HadAM3P from the weather@home distributed climate modelling system 52 give a temperature difference of around 1.5 K between the present and counterfactual (excluding human influence) climate (results not shown here).

A reduction in the HadCRUT4-CW variability is evident prior to 1940 (Fig. 2a) and probably reflects a greater reliance on interpolation over this period. The number of observations in the grid cell increases considerably from 1940, so this is the period regressed against the GWI. Regional anthropogenically forced warming above the 1850-1880 baseline was calculated for 1989-2018 as equalling 95% (85-105%, 5-95% range) of the observed temperature change. This value displays a slight sensitivity to the choice of the present-day period. If a period of ten years is used, following Allen et al. 26 , regional anthropogenically forced warming equals 120% (105-135%, 5-95% range) of the observed value. The regression onto a global index assumes that the proportional contributions of anthropogenic and natural forcing to temperature change are the same at regional and global scales, justified by the absence of independent local factors that could create century-scale local climate trends. Such factors might include volcanic activity with effects localized to the Peruvian climate or substantial local changes in albedo due to land-cover alterations over the past 150 years. No such factors are known in this region 53 .

We generated uncertainty estimates for the anthropogenic contribution to the temperature trends by regressing control simulations of 41 CMIP5 models (with no trends) against the GWI to calculate region-specific scaling factors for GWI. The product of GWI and these scaling factors creates a PDF of the degree that natural variability can project onto the observed temperature trend. The 5-95% range of this PDF is the uncertainty in our estimate of the anthropogenic influence. We note that temperature variability may be higher at the local scale than that simulated in CMIP5. However, in the Cordillera Blanca, the magnitude of interannual variability at grid scale and station data has been found to be similar over a 30 yr period 54 . Further, even given an extreme assumption that local interannual temperature variability is double that of the CMIP5 grid cell, the 5-95% confidence range in the attribution statement broadens only slightly, from 85-105% to 70-120%.

We took annual precipitation values as the accumulation zone mass balance (for 2006-2008) from the nearby Shallap glacier 29 . Previous studies found that most individual meteorological stations do not show significant trends in precipitation over the twentieth century 34 , although some have found an increasing trend in Cordillera Blanca precipitation in recent decades 50 . This glacier is also comparatively insensitive to variations in precipitation (discussed further below) and we did not apply a precipitation trend here.

The HadCRUT4-CW temperature dataset gives a standard deviation of 0.55 K in annual mean temperature for between 1880 and 2018, and the Legates and Wilmot precipitation dataset 55 gives a standard deviation of 0.2 m yr -1 in annual mean precipitation (1900-2010, with no significant trend) 55 . We tested for interannual persistence in the climate data using autoregression modelling 56 . The precipitation data are consistent with no persistence (that is, white noise) but, for the HadCRUT4-CW temperature data, the best-fit autoregressive process is a first-order (that is, AR(1)) process with a lag-1 correlation coefficient of 0.32, which indicates a decorrelation time of around a year. This persistence is probably associated with the influence of El Niño. Our synthetic time series of stochastic climate variability includes this persistence, which enhances the magnitude of natural glacier variability 2 . We also applied the autoregression model to an ensemble of 41 models of CMIP5 control runs and found lesser levels of persistence than observed in almost all models. As interannual persistence enhances the natural variability against which we compare the retreat of Palcaraju glacier, the modelled persistence in CMIP5 suggests our incorporation of interannual persistence is a conservative assumption.

Vertical profile of mass balance. The influence of precipitation on the vertical mass-balance profile is incorporated in our modelling by using the observed mass-balance profile from Shallap glacier (with a similar aspect to that of Palcaraju). Gurgiser et al. 29 ) is approximately zero and the mass balance tapers smoothly to a constant of approximately 1.25 m of water equivalent per year. We join these two domains smoothly with a tanh function of altitude, z:

where z ELA is the ELA and δz is the width of the transition zone (set as 200 m). b(z) is then obtained by integrating equation ( 1) downwards from the top of the domain. Gurgiser et al. 29 Tropical glaciers typically experience precipitation during the ablation season, and frequent changes in snow cover drive strong vertical gradients in albedo, and therefore mass balance 57 . The long ablation season of tropical glaciers also contributes to the strong vertical mass-balance gradient 35,57 , which is nearly constant from year to year. Over two years of observations from the nearby Shallap glacier, Gurgiser et al. 29 reported only a 3% difference in the vertical gradient of the mass balance of the ablation zone, despite a 35% variation in accumulation. Consistent gradients can also be seen, for example, over 10 yr at Uruashraju in Peru 35 , with similar results found for Zongo glacier in Bolivia 58 and elsewhere in the tropics 59 . The near-constant year-on-year shape of the mass-balance profiles means that vertical and horizontal shifts of the mass-balance profile are a useful and efficient approximation to represent climate forcing 30,31,33 .

Climate sensitivity of mass balance. Following refs. 33,35,60 , the impact of temperature anomaly (T′) is implemented by vertically shifting the mass balance by an amount z′ = T′/Γ, where Γ = 5.5 × 10 -3 K km -1 is the observed lapse rate 29,57 . We quantify the impact of a precipitation anomaly by adding a uniform anomaly to the mass-balance profile (a horizontal shift in Fig. 2b). The steep mass-balance gradient in the ablation zone means that mass-balance variability is almost entirely due to temperature fluctuations. The interannual standard deviation of the ELA is 102 m with precipitation and temperature variability included, and 100 m when precipitation variability is omitted. For the 1880 glacier profile (Fig. 3a) this translates to an interannual standard deviation in mass balance of 0.78 m yr -1 and 0.74 m yr -1 , respectively, when precipitation is, and is not, included.

Our calculations, provided in full in the Supplementary Information, show that the ELA rises by 180 m for a 1 K temperature rise, in line with Sagredo et al. 32 (180 m) and Kaser 35 (182 m). By contrast, a 10% increase in annual mean precipitation raises the mass balance by approximately 0.1 m yr -1 , equivalent to a 13 m decrease in ELA. Normalizing the results from other tropical glaciers in the published literature to this same mass-balance perturbation gave ELA decreases of 16 m (ref. 32 ) or 8 m (ref. 35 ). These ELA changes are much lower than those that result from a 1 K increase in temperature. Temperature trends are much more important than the conceivable precipitation trends for the ELA and mass balance of Palcaraju glacier.

We further note that the mass-balance sensitivity to precipitation is reduced by Palcaraju glacier's wet climate, which means that, unlike some tropical glaciers, sublimation contributes only a small portion of the mass budget. Drier glaciers with substantial sublimation tend to have a higher precipitation sensitivity 61,62 . Gurgiser et al. 29 reported that sublimation and evaporation accounts for about 10% of the ablation on Shallap glacier, consistent with the specific energy balance calculations of Rupper and Roe 62 .

Catchment geometry. The geometry of the Palcaraju glacier catchment was obtained from Google Earth (Lake Palcacocha 9° 23′ 42″ S, 77° 22′ 42″ W, elevation 8.8 km). A central-catchment flowline was estimated to obtain a longitudinal profile of the topography. We replaced the lake surface in Google Earth with the observed lake bathymetry 13 , and restored the topography of pre-1941 GLOF moraine at the end of the lake based on the reconstruction in Rivas et al. 13 . The resulting profile is shown in Fig. 3a. The profile does not remove the existing ice in the upper catchment, but as its thickness is only a few tens of metres 63 , it has only a small effect on the profile. Ten approximately evenly spaced swaths along the profile were measured for the catchment width and symmetrized about the central flowline to produce the idealized planform catchment geometry, shown in the insets in Fig. 3a. Numerical model. The numerical model is the same as that used in Roe and O'Neal 6 , based on Oerlemans 64 . The numerical model calculates ice flow from the shallow-ice approximation and Glen's Flow Law, and can be expressed as a pair of equations that are a function of time, t, and distance along the flowline, x:

where h is ice thickness, w is flowline width, F is the downslope flux, b is the mass balance, f d is the deformation constant, ρ is the ice density, g is gravity and dz s / dx is the surface slope. Equations ( 2) and ( 3) constitute a non-linear diffusion equation, wherein the ice flux is a sensitive function of the glacier shape. This numerical model is a useful tool to evaluate the glacier profiles that are consistent with the landscape geometry and the rising ELA associated with the observed climate trend. It also provides a method to estimate the timescale of the glacier response. This response time controls both the magnitude of the natural glacier fluctuations and the glacier's response to a climate trend 5,65 . The flowline model does not incorporate the effect of debris cover, subglacial till or lake-terminus dynamics, which includes ice calving and thermal conduction. Glacier flow may also have altered the subglacial topography which, after the retreat, became the lake bathymetry. As such, the flowline model is not intended or expected to produce an exact simulation of the glacier's evolution.

To estimate the glacier-length response to natural climate variability, we generated 10 kyr of synthetic climate variability with the same statistical characteristics (variance and lag-1-year autocorrelation) as the observations. The resulting glacier-length fluctuations have a standard deviation of 235 m.

Location of the 1880 terminus. In the absence of direct observations, we estimated the late nineteenth century terminus position, L 0 . For our model, L 0 is the also the terminus position in equilibrium with the long-term average pre-industrial climate. Centennial-scale fluctuations about this equilibrium position result from the glacier's response to interannual climate variability (Fig. 4a) 6,65,66 . Our estimate for the magnitude of natural fluctuations is σ L = 235 m. L 0 was necessarily set back from the main terminal moraine as previous work has shown the moraine formed in the seventeenth century 15 , and it has not since been overtopped by natural glacier length fluctuations. Reconstructions from elsewhere in the Peruvian Andes also indicate that glaciers had pulled back from their seventeenth century maxima by the late nineteenth century 42 . It is also unlikely that L 0 was adjacent to the 1940 position: this would require a larger natural variability for the glacier to have extended as far as the seventeenth century moraine, and these larger fluctuations would have taken the glacier further up the valley than the 1940 position, for which there is no evidence here or elsewhere in the region 15,42 . Glaciers with larger natural variability are also more sensitive to climate change 65 , which implies an even more sensitive system than that modelled here. Finally, some glacial retreat in the Cordillera Blanca was observed as temperatures began to rise in the decades prior to 1940 14,[42][43][44] . We therefore set L 0 = 4,250 m. This places Palcaraju glacier's nineteenth century terminus approximately 2σ L from both the seventeenth century moraine and the 1940 position, which is consistent with the glacier fluctuating within this portion of the valley in the eighteenth and nineteenth centuries, but not outside of it (Figs. 3a and4a). Placing L 0 halfway between the seventeenth century moraine and the 1940 moraine is conservative: it allows for the largest magnitude of natural variability and is consistent with the known bounds of the glacier and our direct modelling of σ L . If the 1880 position of the terminus is placed substantially closer to the 1940 front (say, L 0 < 4 km), it also becomes difficult to reconcile that 1880 location with the approximately 40% increase in net ablation rates that accompanies the 0.25 K of anthropogenic contribution to the 1880 to 1940 warming, seen in the insets on Fig. 3.

A comprehensive uncertainty analysis of the change in glacier length that can occur from natural variability alone (ΔL| null ) is possible. Roe et al. 2 used the property of the glaciers' amplification of the signal-to-noise ratio of climate change, and that the amplification is greatest at intermediate glacier response times, to show that:

where γ(t 0 ,τ) is an amplification factor that depends only on the duration of the change (t 0 ) and the glacier response time (τ), and s b | null is the signal-to-noise ratio of the mass balance (b) due to the observed interannual variability (that is, Δb/σ b ). The glacier response time equals the characteristic ice thickness (H) in the ablation zone divided by the (negative) net mass balance at the terminus (b t ) and affects both the amplitude and rate of change of length possible over a given time period, t 0 . For the glacier's 1880 geometry, we estimated the response time as τ ¼ H �bt  150 m=3 m yr �1  50 yr I . We have model estimates of τ and σ L , but they are uncertain. We represent uncertainty in τ with a gamma function with 95% confidence bounds of 25-75 yr (Fig. 4d). Similarly, for σ L , we give 95% confidence bounds of 115-345 m (Supplementary Fig. 2c). s b | null is taken from the observed interannual variability and follows the Student's t-distribution 2 . Through equation ( 4), these PDFs yield PDFs for trends of length 60 and 140 yr in the absence of a climate change (Fig. 4e,f).

Equation ( 4) is derived in Roe et al. 2 based on a three-stage model of glacier dynamics, whose parameters are governed by the glacier's geometry and climatic setting 65 . The model has been shown to accurately emulate glacier fluctuations in numerical models of ice-flow dynamics [65][66][67] . The three-stage model is linear, and although we do not expect linearity to apply over the whole valley profile, our goal was to characterize the much smaller fluctuations due to natural variability. The synthetic mass balance time series we applied to the numerical model (Fig. 4a) is normally distributed, so the resulting glacier fluctuations are also approximately normally distributed, which supports the linearity approximation. GLOF hazard assessment. We used two methods to identify dangerous glacial lakes. The method proposed by Wang et al. 46 was developed to assess settings in which the likely GLOF mechanism would be a moraine-overtopping wave induced by an ice avalanche, as is the case with Lake Palcacocha. The GLOF hazard is quantified using the glacier area, the distance between the lake and glacier terminus, the slope between the lake and glacier, the mean slope of the moraine dam and the glacier-snout steepness. As Lake Palcacocha is in direct contact with Palcaraju glacier, we followed Wang et al. 46 and used the gradient of the glacier snout, defined as the glacier area within a 500 m horizontal distance of the terminus, in place of the slope between the lake and glacier. The area of the glacier is taken from the Global Land Ice Measurements from Space database 68 (the value is for 2015). The mean slope of the moraine dam before and after the 1941 GLOF (for the pre-industrial and present-day calculations, respectively) was taken from Somos-Valenzuela et al. 12 , and from the slope of the glacier snout from our model reconstructions of the 1880 and present-day glacier geometries. Our assessment gives the same numerical value for the present GLOF hazard as that of Emmer and Vilímek 11 , using this methodology.

We supplemented these calculations with an independent assessment adapted from Bolch et al. 47 . This methodology uses 11 weighted stability parameters, which relate to characteristics of the lake, lake surroundings and adjacent glaciers, and the potential impact on downstream areas. The parameters were ranked according to their influence on GLOF occurrence. We excluded two parameters from our analysis (glacier terminus slope and glacier velocity) as our analysis concerns the present-day GLOF hazard and these factors are used to identify glacial lakes for which a future rapid expansion may occur, which increases the hazard. We recalculated the weight factors for the remaining nine parameters according to the method provided in Bolch et al. 47 (the values are provided in Supplementary Table 3), which also instructs that flash-flood risk is set to zero if a debris flow is expected. We estimated the nineteenth century lake area as approximately 150,000 m 2 based on our modelling estimates of glacial retreat prior to 1941 and the estimated area of 303,000 m 2 immediately prior to the 1941 GLOF 69 . The lake area was 518,000 m 2 in 2016 23 . This method assesses hazard as 'high' if the four most important factors apply or if a combination of factors that, combined, is equal to the sum of the weights of the four most important factors apply.