Little Ice Age climate reconstruction from ensemble reanalysis of Alpine glacier fluctuations

Mountain glaciers sample a combination of cli- mate fields - temperature, precipitation and radiation - by accumulation and melting of ice. Flow dynamics acts as a transfer function that maps volume changes to a length re- sponse of the glacier terminus. Long histories of terminus po- sitions have been assembled for several glaciers in the Alps. Here I analyze terminus position histories from an ensemble of seven glaciers in the Alps with a macroscopic model of glacier dynamics to derive a history of glacier equilibrium line altitude (ELA) for the time span 400-2010 C.E. The re- sulting climatic reconstruction depends only on records of glacier variations. The reconstructed ELA history is similar to recent reconstructions of Alpine summer temperature and Atlantic Multidecadal Oscillation (AMO) index, but bears little resemblance to reconstructed precipitation variations. Most reconstructed low-ELA periods coincide with large ex- plosive volcano eruptions, hinting at a direct effect of vol- canic radiative cooling on mass balance. The glacier ad- vances during the LIA, and the retreat after 1860, can thus be mainly attributed to temperature and volcanic radiative cool- ing.


Introduction
Proxy-based climate reconstructions are an important source of information on past climate variability and help to assess whether the climate change presently observed is unprecedented in a long-term context (IPCC, 2007). Some of these proxies, notably tree-ring data, give a detailed picture of summer temperature variations (e.g., Mann et al., 2009;Büntgen et al., 2006Büntgen et al., , 2011, while historic records provide more qualitative, but very detailed information (e.g., Pfister, 1999;Casty et al., 2005). Records of glacier length changes (GLCs) provide supplementary independent information but are caused by fundamentally different processes.
GLCs have been used in a qualitative manner to confirm climate reconstructions from other proxies (Guiot et al., 2010) or compare them with solar forcing reconstructions (e.g., Beer et al., 2009;Nussbaumer et al., 2011). However, such comparisons suffer from the delayed response of glacier length to climatic or radiative forcing, and are only valid on long timescales.
Mountain glaciers adjust their length to a changing climate at a rate which depends on bedrock geometry, ice thickness distribution and local mass balance distribution. Glaciers can therefore be regarded as climate sensors with geometrydependent susceptibility and response time. This relationship has been exploited to infer past climate from GLCs with the help of glacier models of different complexity. Analytical approximations (Callendar, 1950;Nye, 1963Nye, , 1965Klok and Oerlemans, 2003;Roe and O'Neal, 2009), neural networks  and numerical ice flow models (Greuell, 1992;Schmeits and Oerlemans, 1997) have been applied to either single glaciers, an ensemble of glaciers in the same climate region (Oerlemans, 2007), or glaciers distributed around the globe (Hoelzle et al., 2003;Oerlemans, 2005;Leclercq et al., 2011;Leclercq and Oerlemans, 2012).
Several studies showed that numerical models of reduced complexity give volume and length changes that are similar to those calculated with models solving the full equations (e.g., Leysinger-Vieli and Gudmundsson, 2004;Lüthi, 2009;Oerlemans, 2011). A simple model of a mountain glacier describing the delayed response of GLC to volume change is  Zumbühl and Holzhauser (1988) Unterer Grindelwald 393-1983Messerli et al. (1975; Zumbühl (1980) Holzhauser andZumbühl (2003) a dynamical system in the variables "length" L and "volume" V derived in Lüthi (2009). Essentially, the glacier behaves as a slightly over-damped harmonic oscillator in one of the macroscopic variables L or V around a steady state (Harrison et al., 2003). The advantage of using a macroscopic model is that no knowledge of the usually complex and insufficiently known glacier bed geometry or basal sliding is required. Instead, two model parameters are determined with an optimization procedure, which yields the LV -model with the most similar length response under a given ELA history. Such a dynamically equivalent simple model (DESM) is useful to systematically analyze GLC records , and to infer glacier volume change from GLCs .
Here I use seven very long glacier terminus position records from the Alps to derive glacier equilibrium line altitude (ELA) for the period 400-2010 C.E. The ELA history causing the observed GLCs is determined with an optimization procedure which employs the macroscopic LV -model for each of the seven glaciers. The resulting glacier ELA variations are completely independent of any other proxy data, and are indicators of summer temperature, winter precipitation and radiation (Vincent et al., 2004;Ohmura et al., 2007;Huss et al., 2009;Thibert et al., 2013). The reconstructed ELA histories are then compared to other proxy-based reconstructions of European summer temperature, total solar irradiance, and volcanic cooling.

Length change data
In the Alps, GLCs have been systematically monitored since about 1860 when most glaciers receded from their last maximum extent of the Little Ice Age cold period (LIA, ca. 1550(LIA, ca. -1850. Older reconstructions are based on photographs, paintings, written records, and the positions of trees destroyed by glacier advance, and dated by dendrochronology or the carbon-14 method. The longest such reconstruction has been assembled for Grosser Aletsch glacier, and similar records exist for the glaciers Gorner, Rhone, Rosen-laui, Unterer Grindelwald, Mer de Glace and Bossons (data sources are listed in Table 1). These historical GLC series have been complemented with recent data from the Swiss Glacier Monitoring Network (Glaciological Reports, 1881 and other sources (Vincent et al., 2005, A. Bauder, personal communication, 2012.

Macroscopic glacier model
Climatic interpretation of GLC has to take into account the combination of direct and delayed response of the glacier terminus to climate forcing. The direct response is caused by the difference between the rates of ice melt and mass transport to the terminus. Ice melt varies rapidly with changes of temperature and radiation, while mass transport varies slowly due to upstream geometry changes, which affect ice flow patterns. A glacier model represents the transfer function from volume change to length change, and depends on glacier size, bedrock geometry, and surface mass balance distribution.
The model used here represents glacier dynamics on a macroscopic level, and results in a dynamical system in the variables "length" L and "volume" V for unit glacier width, referred to as the LV -model (Eq. 40 in Lüthi, 2009): Here,ġ = @ḃ @z denotes the vertical gradient of local mass balance rateḃ (in units of velocity, i.e., meter ice equivalent per year), and s = tan is bed slope. Equation (1b) is a relaxation equation for the current glacier length L with time constant ⌧ a (in years; ⌧ a should not be confounded with the volume timescale ⌧ v ). The steady-state length for the current volume V is determined by a volume-length scaling relation of the form V = a L µ with parameters µ = 7/5, and a which depends onġ and s (Lüthi, 2009, Eq. 21), and is continually updated for the current glacier length (Lüthi, 2009, Eqs. 37 and 38). The dynamical system (Eq. 1) is driven by a forcing in the term Z(t) = z 0 z ELA (t), where z 0 is the highest point  (Table 1 for data sources). of the bed and z ELA (t) is the time dependent elevation of the equilibrium line.

ELA reconstruction
The LV -model (Eq. 1) was solved forward in time with an ODE solver (LSODA from ODEPACK, Hindmarsh, 1983), and evaluated on a yearly time step. To drive the LVmodel and hence calculate glacier response, an ELA history z ELA (t) is prescribed as a forcing function, changing the Z(t) = z 0 z ELA (t) term of Eq. (1a). In the LV -model, local mass balance rate is given as linear function of elevatioṅ b(z) =ġ(z z ELA ) with a constant vertical gradient of local mass balance rateġ. A value ofġ = 0.008 yr 1 was assumed for this study, which is a good average for the central Alps (e.g. Chen and Funk, 1990;Huss et al., 2008; Table 3); usinġ g = 0.007 yr 1 changed the results only marginally.
The climate history is represented by constant ELA values within fixed size intervals (10 years between 0 and 1200 C.E., 5 years between 1200 and 2010 C.E.). These 283 ELA values are unknowns to be determined in an optimization procedure. Additionally, two parameters (slope s and vertical extent Z 0 = Z(t = 0)) are sought for each of the seven LVmodels representing the seven glaciers. Optimization was performed with two algorithms, leading to very similar re-sults: a reversible jump Markov chain-Monte Carlo (RJ-MCMC) algorithm, and the BOBYQA optimizer (Powell, 2009;Johnson, 2012). The cost function for the optimization is constructed by summing the squares of differences between modeled and measured GLCs. A large penalty was added to the cost in case of violation of one-sided constraints on glacier length (e.g., living trees or terminus moraines), and no cost was added if no constraints are available. The model runs were initialized with random glacier parameters, and zero (or random) ELA variations.

Reconstructed length changes
Modeled terminus positions calculated with the bestperforming ELA histories are compared to the historical records in Fig. 1 (Table 1 for   individual terminus positions since 1600 C.E. when dense data are available. Root mean square (rms) differences are smaller than 20 m. Most absolute differences are smaller than 150 m, but are occasionally as high as 300 m. Since the LVmodel geometry with a single slope is highly idealized, perfect agreement between observed and modeled GLCs cannot be expected. Terminus positions are strongly influenced by details of terminus bedrock geometry such as cliffs (glaciers Mer de Glace, Unterer Grindelwald and Rosenlaui) or terminus areas with strong topographic breaks (Rhone). By using an ensemble of glaciers with different geometries and response characteristics, the dependence of the overall result on such geometric peculiarities is reduced. Table 2 lists the optimized LV -model parameters of the best-fitting DESMs, and their volume and length timescales. Vertical elevation differences, Z 0 , and bed slopes, s, differ considerably from the values of real glaciers. As discussed in  these quantities include the effects of basal motion and terminus geometry, which are ignored in the simplistic LV -model dynamics and geometry.
The modeled GLCs obtained with the reconstructed ELA histories are robust against omission of length constraints and variations of glacier parameters (changing s by 5 %, and Z 0 orġ by 10 % leads to 1 % changes in length change variability). During several time periods with small glacier extent no terminus position data are available and modeled terminus positions assume a wide spread of values (around 700, 900, 1000 and 1400-1530 C.E., hidden red lines in Fig. 1). Consequently, no information on ELA variations during these periods can be extracted, and the modeled terminus positions and reconstructed ELAs can assume a wide variety of values within these periods.
The low reconstructed ELAs during the period 1300-1350, only documented from trees killed by the advancing Gorner glacier, lead to large modeled advances of several glaciers that exceed the 1850 terminus positions. Indeed, on the Rhone glacier several moraine ridges from the early medieval period are documented, which could correspond to this reconstructed 1350 maximum position (Zumbühl and  Holzhauser, 1988). The slowly reacting glaciers Gorner and Grosser Aletsch on the other hand only reach terminus positions that are less advanced than their 1850-1860 maximum.

Reconstructed ELA history
The reconstructed ELA histories leading to the best agreement between modeled and measured GLCs are shown in Figs Periods of high ELA (high summer temperature and radiation, low solid precipitation), corresponding to decreasing or small glacier extents, are mostly missing in the reconstruction before 1600. This is due to the lack of data on shorter glacier extents, traces of which might still be buried under the glaciers. The period with highest ELAs begins after 1980, exceeding the reconstructed values of the preceding 1600 years. Since present-day glacier geometries have not yet fully adjusted to the climate of the last decades (depending on glacier response time) the reconstructed ELA during the last decades might be considerably higher. The reason for this incomplete recent reconstruction is the delayed reaction of the terminus to mass balance changes.

Comparison with instrumental data
The reconstructed ELA history is compared to instrumental temperature data in Fig. 2. The agreement of summer (JJAS) temperature with the high-Alpine Grand Saint Bernard meteorological station (2500 m a.s.l., Begert et al., 2005) and with the HISTALP CRSM Alpin multi-station reconstruction (Auer et al., 2007) in Fig. 2b is reassuring. To obtain similar variability of the records, the sensitivity factor between reconstructed ELA and temperature is f s = 80-100 m K 1 (i.e., 80-100 m ELA change per Kelvin temperature change). Similar values are reported in the Fig. 2. ELA variation reconstruction leading to best agreement with length change records from seven glaciers (Fig. 1). Modeled ELAs are shown with red and blue colored areas. Instrumental records are shown for comparison in panels (B) and (C) as 5-year running means of summer (JJAS) temperature (data from MeteoSwiss (Begert et al., 2005), Casty et al. (2005) and HISTALP (Auer et al., 2007)). Panel (D) displays winter and spring precipitation as 5-year running means (Pauling et al., 2006). (A) and vertical stripes show timing and magnitude of the Northern Hemisphere stratospheric sulfate loading from explosive volcanism (Gao et al., 2008;Arfeuille et al., 2014). literature: 97 m K 1 (Thibert et al., 2013), 113 m K 1 (Oerlemans, 2010), and 60-70 m K 1 (Vincent, 2002).
Comparisons with the HISTALP CRSM NW (Auer et al., 2007) and the Casty et al. (2005) reconstruction of European temperature, reaching back to 1768 and 1500 C.E., are shown in Fig. 2c. The ELA reconstruction shows similar temperature variability after 1820, but is considerably lower before (by 100 m, or 1 K), while the timing of high and low temperature events tends to agree. Such a discrepancy might be caused by a suspected warm bias in the early instrumental data (Frank et al., 2007;Böhm et al., 2010).
To evaluate the role of precipitation, a comparison with the gridded reconstruction form Pauling et al. (2006) is shown in Fig. 2d. Average precipitation values of winter and spring precipitation from 12 grid cells over the Alpine area are displayed. The low-ELA phases of 1712-1717, 1747-1752, 1912-1922 and 1952-1957 coincide with high winter precipitation, while the periods 1712-1717, 1767-1777 and 1952-1957 coincide with high spring precipitation. None of the other precipitation variations are reflected in the reconstructed ELAs. Summer precipitation is mostly uncorrelated with the reconstructed ELAs and therefore not shown.

Comparison with proxy records
The reconstructed ELA history was compared to many long climate reconstructions of mid-European temperature from proxy data such as tree-rings, lake sediments, ice cores and speleothems. Since glacier mass balance is controlled by a combination of processes, one cannot expect one single proxy record to agree with the ELA history. In addition, climate variability at high elevation might differ significantly from the low elevations from which most proxy records are recovered. Figures 3 and 4 compare the reconstructed ELA history to a selection of proxy records that show, at least during limited time spans, close similarities. Reconstructions of Atlantic Multidecadal Oscillation (AMO, Fig. 3c; Mann et al., 2009) and a Northern Hemisphere temperature history (annual mean, multi-proxy; Fig. 3d; Mann et al., 2009) show similar long-term behavior, variability and individual cold events after 1200, but are considerably warmer before. The sensitivity factor from ELA to temperature variation is twice that of the instrumental data (Table 3), presumably due to hemispheric averaging of the records. A recent multi-proxy annual-mean temperature history for Europe (PAGES 2k consortium, 2013, Fig. 4e;) is very similar to the ELA reconstruction throughout the considered 1600 years. It shows most of the early cold events (550, 800, 1100 C.E.) that are absent in other proxy records, but for which marked glacier advances are documented (Fig. 1). Tree-ring-based temperature histories of the Alps ( Fig. 4f; Büntgen et al., 2006) and central Europe ( Fig. 4g; Büntgen et al., 2011) show an overall close agreement with the ELA reconstruction. Tree growth is predominantly a proxy for summer (JJAS) temperature, which coincidentally is the same period which affects glacier melt, and thus glacier mass balance variability (e.g., Huss et al., 2008). Figure 3b shows that reconstructed ELA variations resemble total solar irradiance (TSI; Steinhilber et al., 2009) during the period 1700-1950, but bear little similarity before 1700. The cooling effect of explosive volcanism can be evaluated with the top panels and gray vertical bars in Figs. 2-4, indicating the Northern Hemisphere stratospheric sulfate load from big volcanic eruptions (Gao et al., 2008;Arfeuille et al., 2014, NH, v2). It is noteworthy that most of the low-ELA phases during the LIA coincide with, or are preceded by, high stratospheric sulfate loads from volcanic eruptions. This aspect is investigated in the discussion below. Figure 1 shows the remarkable result that one single ELA history (Fig. 3) causes long-term GLCs of seven Alpine glaciers with very different geometries that mostly agree with the documented record. This result justifies a posteriori the important assumption that all GLCs are caused by the same ELA history. This assumption is also supported by studies of glacier mass balance variability (Vincent et al., 2004;Huss et al., 2009). Even if all glaciers are in the same mountain range and within 130 km of each other, local climate, and especially precipitation, vary on short spatial scales, while their variability is similar (e.g., Frei and Schär, 1998;Casty et al., 2005). Long-term instrumental records show that variations in temperature are closely linked over the Alps (e.g., Casty et al., 2005;Auer et al., 2007).

Discussion
Earlier climate reconstructions from GLCs (Klok and Oerlemans, 2003;Oerlemans, 2005;Leclercq et al., 2011;Leclercq and Oerlemans, 2012) used a linear response model to derive temperature for each glacier individually from smoothed length data. By contrast, the presented method uses a more accurate representation of macroscopic glacier dynamics (LV -model, Eq. 1), unaltered and much longer GLC records, and optimizes on a single ELA history constrained  (Gao et al., 2008). (E) multi-proxy temperature reconstruction (PAGES 2k consortium, 2013). (F, G) temperature reconstructions from tree rings (Büntgen et al., 2006(Büntgen et al., , 2011. All temperature curves are smoothed with a 7-year running mean.
by the responses of seven glaciers with very different geometries. An important aspect of the presented method is that terminus positions documented for some glaciers can be used to infer missing data of the remaining glaciers. A sufficiently long calibration time span with well-documented glacier advances and retreats is a prerequisite to obtaining good parameters for the DESM glacier models. Obviously, modeled GLCs still differ from reality by the simplistic model dynamics and the assumed geometry which cannot represent the dynamical effects of small-scale terminus topography. However, the overall timing and range of GLCs are stable.

Warm periods
The ELA values during warm periods with significant glacier retreat cannot be inferred from the terminus position record because of lack of data on minimum extents for the considered glaciers (which might still be buried under ice). Extended warm periods with glacier extents equal to, or smaller than today are possible within the gaps in the record, i.e., between 650-800, 850-950, 980-1050, 1200-1280, and 1400-1550 C.E. Recent findings of washed-out peat and wood from pro-glacial forefields of other Alpine glaciers indicate glacier extents smaller than the present between 690-780 (Hormes et al., 2001(Hormes et al., , 2006Joerin et al., 2006). No similar findings are reported for the other periods.

Climatic interpretation
For the climatic interpretation of the reconstructed ELA history it is important to discern the different climate components. Summer mass balance, the balance component with high interannual and interdecadal variability (Huss et al., 2008), depends mainly on temperature and to a lesser degree on radiation (as measured at remote stations; Ohmura et al., 2007), while winter mass balance is dominated by solid precipitation. The influence of precipitation on the ELA is relatively small: 1 K temperature change is equivalent to a precipitation change of 350 mm yr 1 w.e., or a radiation change of 7 W m 2 (Ohmura et al., 1992). The importance of atmospheric radiation variations (solar dimming through aerosols; changes in cloudiness) on mass balance has been highlighted by Huss et al. (2009). In a longer-term context, summer temperature and precipitation have been identified as dominant controls on Aletsch glacier mass balance (Steiner et al., , 2008 Several recent studies state that GLCs are closely related to total solar irradiance (TSI), although the similarity is limited to several time spans (e.g., Beer et al., 2009;Nussbaumer et al., 2011). Since GLCs are mostly a delayed reaction to climate forcing, a more direct relationship should exist between ELA variations and shortwave radiation. The comparison with TSI in Fig. 3b does not support the hypothesis of a direct influence of solar variability on glacier mass balance, except for the limited time span 1700-1950, where the similarity is striking. After 1950 the TSI declines, whereas ELA strongly rises. Before 1700 the two quantities seem mostly unrelated. The above reasoning does not exclude the potentially important radiation modulation by atmospheric processes (aerosols, moisture, cloudiness) which is likely important for glacier mass balance (Ohmura et al., 2007;Huss et al., 2009). Since no long-term reconstructions of these parameters exist, their influence cannot be meaningfully investigated with the presented ELA reconstruction.

Volcanic cooling
A potentially important cause for the episodic and rapid glacier advances, reflected in the reconstructed low ELA periods, is radiative summer cooling by big volcano eruptions. Stratospheric sulfate aerosols from large explosive tropical volcano eruptions cause northern hemispheric summer cooling and altered circulation patterns (e.g., Crowley, 2000;Robock, 2000;Zanchettin et al., 2013). They cause incoming solar radiation reduction by several W m 2 for 2 years (e.g., Robock, 2000;Weber, 2005;Fischer et al., 2007;Zanchettin et al., 2012), and an effect of 0.2 K cooling lasting about 4 years in central Europe (Esper et al., 2013b). The full range of cold periods caused by volcanic events might not even be recorded in tree rings (Mann et al., 2012). Increased summer precipitation in the following cool summers was attributed to the eruptions (Wegmann et al., 2014) which likely leads to considerably more positive mass balances. Longer-term persistence of cold summers is likely through sea ice/ocean feedback mechanisms, possibly even explaining the onset of the LIA cold period (Miller et al., 2012;Schleussner and Feulner, 2013) although this is a topic under debate (Esper et al., 2013b).
Northern Hemisphere stratospheric sulfate loads reconstructed from ice core data (Gao et al., 2008) are shown in panels A of Figs. 2-4. It is striking that most reconstructed low-ELA phases during the LIA coincide with, or are closely preceded by, high stratospheric sulfate loads (marked with gray vertical bars). Thus volcanic cooling might explain most of the short-lived, rapid LIA glacier advances of the glaciers Bossons, Mer de Glace, and Unterer Grindelwald (Fig. 1). Low-ELA phases since 1580 (in 5-year steps from the reconstruction) are compared to eruptions in Table 4 and Fig. 2. All Table 4. Phases of very low reconstructed ELA between 1580 and 1920, and explosive tropical volcano eruptions causing significant aerosol loading in the Northern Hemisphere (NH). The 1t row indicates the number of years of delay (+) or advance ( ) with respect to the eruptions. The estimated NH stratospheric aerosol loading (based on ice cores, Gao et al., 2008) is given in the last column (data from Table 1 in Arfeuille et al., 2014).  2 19122 19222 +0 1912 of these cold phases are within ±4 years of major eruptions, with the exception of the very large 1783 Laki event, which precedes reconstructed low ELAs by 6 years. This eruption was mainly tropospheric and caused a very cold year 1783/84 in Europe (e.g., Robock, 2000). The hypothesis of a volcanic cooling effect on glacier mass balance is supported by Table 4, but the occurring time shifts of several years have to be explained. While the accuracy of the timing of big volcano eruptions is good (layer counting in ice cores), reconstructed glacier terminus positions are less well constrained. Old paintings are sometimes not accurately dated, and radiocarbon ages of trees have sometimes large error bars. Moreover, tree trunks are recovered from a variety of positions with respect the glacier snout, which itself is not of constant shape, so that another source of dating error results. Furthermore, the ELA reconstruction uses bins of 5 years with constant ELA, such that the dates of reconstructed ELAs are of the same order at best. Whether such effects explain the time lags between volcanoes and cold ELA phases could be assessed by forward modeling of GLCs with a combination of climate time series as forcing functions. Such work is currently in progress, and hints at an important role of volcanic forcing.
Before 1500 C.E. no close correlation between volcanic activity and reconstructed ELAs, respectively GLCs, is discernible in Fig. 3. Whether this is due to the much sparser glacier length data, only recording large advances and thus leading to a less-constrained ELA history, or to other causes, cannot be assessed. Among the low-ELA phases without accompanying volcanic event, those of 510-580, 810-830, and 1330-1340 all correspond to cold phases in the temperature reconstructions from Büntgen et al. (2011) andPAGES 2k consortium (2013).
The absence of any recorded glacier advance after the very large 1258 volcanic eruption (Oppenheimer, 2003) is noteworthy. Although the dimming effect of volcanic aerosols only pertains to a few years, advances of several hundred meters of the fast-reacting glaciers (Bossons, Mer de Glace, Rosenlaui) were recorded following each volcanic event after 1600. The lack of traces of a glacier advance after 1258 hints at small glacier extents during the period 1200-1300. Huss et al. (2010) concluded from a study of 20th century glacier mass balance that the Atlantic Multidecadal Oscillation (AMO) strongly influences glacier mass balance in the Alps. Figure 3c shows the comparison between a long AMO series (Mann et al., 2009) and the reconstructed ELA history. After 1250 C.E. the timing of warm and cold phases and the variability is very similar, although the low ELA events still would have to be explained by volcanic cooling. Our ELA reconstruction could thus be interpreted to support the hypothesis that the AMO influences mass balance of Alpine glaciers.

The riddle of the Little Ice Age glacier advance
Several studies concluded that the rapid glacier advances during the Little Ice Age (LIA, ca. 1550-1850 C.E.), and the last maximum extent around 1850-1865 cannot be explained using homogenized instrumental temperature data (e.g., Nesje and Dahl, 2003;Vincent et al., 2005;Nussbaumer et al., 2011). Likely processes leading to this rapid glacier growth are increased accumulation due to 25 % higher winter precipitation than present (Vincent et al., 2005), or reduced melt due to lower incoming shortwave radiation (Miller et al., 2012). More qualitatively, Nussbaumer et al. (2011) conclude that "a combination of low solar forcing, frequent and strong volcanic eruptions, and dynamic effects due to internal variability of the climate system led to the prominent glacier advances during the LIA". The rapid recession of the glaciers after 1865 has been explained by glacier surface darkening, and thus albedo reduction, by increasing release of industrial black carbon (Painter et al., 2013).
The agreement between our ELA reconstruction and several reconstructions of temperature (Fig. 4e, PAGES 2k consortium, 2013; f, Büntgen et al., 2006;g, Büntgen et al., 2011) and AMO index (Fig. 3c, Mann et al., 2009) show that additional large changes in winter precipitation are not needed to explain the rapid glacier advances, although higher summer precipitation has been inferred between 1700 and 1810 from tree-ring reconstruction (Büntgen et al., 2011), and wetter winters were frequent between 1800-1850 (Paul-ing and Paeth, 2007). The frequent short phases of very low ELA during the LIA can be mostly explained by the effect of volcanic eruptions on incoming radiation. Also, the rapid melt-back of the glaciers after 1860 can be explained by temperature alone, without additional albedo changes due to black carbon (Painter et al., 2013).
The question remains why before 1850 the homogenized instrumental temperature records (Casty et al., 2005;Auer et al., 2007) are considerably (⇠ 0.5-1 K) higher than the recent temperature reconstructions of Buentgen (tree rings) and PAGES (multi-proxy). One explanation is that early instrumental measurements were biased because of inadequate shading from solar radiation of the thermometers (Frank et al., 2007;Böhm et al., 2010). Another possible explanation is the lack of instrumental temperature data from the Alps before 1850, which requires the extrapolation of data from distant stations.

Conclusions
The presented method to reconstruct ELA variations from a set of GLC records yielded a new 1600-year history, which is independent of other instrumental or proxy data. Through the simultaneous optimization on the GLCs of seven glaciers with one ELA history, the result is less susceptible to geometry-dependent peculiarities of terminus dynamics of individual glaciers. Resulting GLCs automatically fill gaps in individual length change records, and recorded variability from fast-reacting glaciers is transferred to slowly reacting glaciers.
The obtained ELA history agrees well with tree-ringderived-and multi-proxy temperature reconstructions. Also, a close correspondence with the AMO index between 1250 and 2010 could be observed. Long-term and short-term variability of these reconstructions are reproduced with the glacier record. This shows that analysis of temperature records with independent GLC records is a useful tool to independently validate these proxy-based reconstructions.
Additionally, the reconstructed ELA history shows cold phases only partially recorded in the proxy records, most of which closely correspond to explosive tropical volcanic eruptions during the LIA period. This hints at a strong direct influence of stratospheric radiation scattering from volcanic aerosols on the energy balance on the glacier surface. The dependence of summer mass balance on radiative forcing and the cooling due to changing stratospheric volcanic aerosols are not yet well understood.
The documented rapid advances of Alpine glaciers during the LIA, and their rapid retreat after 1865, can be reproduced with ELA changes that correspond to temperature reconstructions, if taking into account volcanic cooling and increased summer precipitation after large eruptions. Alternative explanations stressing the importance of increased winter precipitation or increased ablation area albedo changes www.the-cryosphere.net/8/639/2014/ The Cryosphere, 8, 639-650, 2014 due to industrial black carbon, are not needed to explain the observations. Work to identify the relative importance of different climate fields for the resulting GLCs with help of macroscopic glacier models is currently under way. An additional avenue for further analysis are the glacier volume changes that the LV -model calculates simultaneously. After calibration with measured volume change data, such model results could be used to quantify long-term runoff changes from glaciers. Such an analysis might decide a 150year-old dispute as to whether glacier melt is to blame for frequent flooding and very high lake levels of Lac Léman (Lake Geneva) in the 1870-1880s, which was the main reason for the long and detailed investigation of glacier changes.