Modelled glacier equilibrium line altitudes during the mid-Holocene in the southern mid-latitudes

Introduction Conclusions References Tables Figures


Introduction
Deciphering the climate signals and glacial history of the mid-latitudes of the Southern Hemisphere (Fig. 1) during the Holocene is key to unravelling the mechanism of climate change that occurred during this period.During the last ∼ 11 500 years, a series of intervals of rapid climate changes occurred worldwide (Mayewski et al., 2004).Reduction in temperature and/or increases in precipitation during these periods have been recorded as multiple glacial advances in different areas of the planet.Solomina et al. (2015) recently provided a global review of Holocene glacier activity.Periods of renewed glacial activity, known as neoglaciations (Porter and Denton, 1967), were initially identified in the North- ern Hemisphere.However, during the last decades, numerous studies have shown evidence of glacial advances, as well as climate variability during this period in the Southern Hemisphere (approximately between 35 and 55 • S).Most of these studies have focused on Patagonia (e.g.Clapperton and Sugden, 1988;Porter, 2000;Rodbell et al., 2009;Strelin et al., 2014) and in the Southern Alps in New Zealand (e.g.Gellatly et al., 1988;Porter, 2000;Schaefer et al., 2009;Putnam et al., 2012;Kaplan et al., 2013).
In the South Island of New Zealand (Fig. 1), on the other hand, notable periods of glacier still stand or re-advance occurred during the early to mid-Holocene, as well as during the last millennium (Schaefer et al., 2009;Putnam et al., 2012).It appears that these centennial-scale glacial events have been superimposed on a long-term trend of decreasing ice extent that persisted for the entire Holocene (Schaefer et al., 2009;Putnam et al., 2012;Kaplan et al., 2013).
In this context, it is clear that several aspects of the neoglacial chronology of the southern mid-latitudes (35-55 • S) are still inadequately understood, and more de-tailed chronologies are needed.Particularly relevant for this study is the lack of agreement regarding the timing of the neoglaciations in the southern mid-latitudes (Porter, 2000).
Understanding the climate and glacial history of the southern mid-latitudes is a prerequisite for testing hypotheses regarding the origin and propagation of palaeoclimate signals, the coupling of the ocean-atmosphere in the extra-tropics, and the interaction of low-and high-latitude climate controls on hemispheric and global climate (Fletcher and Moreno, 2012;Moreno et al., 2010;Rojas et al., 2009;Putnam et al., 2012).The mid-Holocene represents a key moment in our late climate history.This period, within the current interglacial cycle, had important differences in orbital parameters with respect to the present conditions and was devoid of influences from late-glacial climate change (Braconnot et al., 2007).Although recent work has demonstrated that orbital forcing may not have played a critical role in glacier behaviour during cold phases of the last glacial cycle (Doughty et al., 2015), the climatic boundary conditions at that time were very different than during the Holocene.Considering the uncertainty in the timing of the beginning of the neoglaciation, but that geologic evidence suggests that glaciers were larger than present between 8000 and 6000 years BP in both regions (Kaplan et al., 2013;Schaefer et al., 2009;Putnam et al., 2012;Douglass et al., 2005;Harrison et al., 2012), the aim of this work is to undertake a firstorder modelling study to assess whether the climatic conditions during the mid-Holocene would permit the existence of more extensive glaciers than during the pre-industrial (PI) period.The PI period was chosen as the reference period instead of 20th century climate because of the availability of modelling output to drive the mass balance model.
At the mid-Holocene (MH hereafter), orbital differences resulted in negative insolation anomalies from November through March and positive anomalies from June to October in southern mid-latitudes (Rojas and Moreno, 2011).We expect that these orbital/insolation differences had a major impact on the glacial extent and especially in the equilibrium line altitudes (ELAs) of glaciers.The ELA is a climatesensitive parameter marking the location on a glacier where accumulation of snow is exactly balanced by ablation (net surface mass balance equals zero; Porter, 1975;Cogley et al., 2011).ELA can be estimated by fitting a curve to data representing surface mass balance as a function of altitude or a mass balance profile (Cogley et al., 2011).It is also possible to estimate a climatic ELA as an average of ELA during 30 or more years (Bakke and Nesje, 2011).In this sense the most appropriate definition for a climatic ELA is the steadystate ELA or long-term ELA that could be estimated by modelling (Cogley et al., 2011).Fluctuations of ELA have been extensively used in palaeoclimatic reconstructions because the ELA is primarily controlled by temperature and precipitation (Porter, 1975;Sagredo et al., 2014).
In this paper, we explore the differences in the climatic ELA between MH and PI conditions, using a degree day model with data based on Paleoclimate Modelling Intercomparison Project Phase 2 (PMIP2) climate models output in the Southern Alps in New Zealand and the Patagonian Andes in South America (Fig. 1).Therefore this study goes a step further towards linking PMIP2 model simulations, and hence orbital forcing, with the MH glacial record by explicitly calculating the regional ELAs in Patagonia and New Zealand at 6000 years BP.The ELA difference between the PI and MH provides information on the scale of glacier change at this key time.

Experimental setup
With the aim of obtaining comparable results we use a glacier mass balance model forced with PMIP2 models for both periods, MH and PI.In the next sections we explain both the model and the data.The general approach consists in resizing all PMIP2 models' outputs to a resolution of 0.5 • using linear interpolation.Due to the coarse resolution of the PMIP2 models, and the regional nature of this study, we used the ELA as a general indicator of glacier behaviour as we are not considering individual glaciers and their specific responses to climatic variations.For each grid point we obtained surface mass balance as a function of altitude.From this mass balance profile we thus obtained a regional ELA.
We applied the same procedure for both time slices (mid-Holocene and pre-industrial) and focus in the difference between the two periods.

Glacier mass balance model
We applied a simple glacier mass balance model to explore the regional differences in the ELA between the MH and the PI in the southern mid-latitudes.
Degree day models, as used in this study, have previously been applied for palaeoclimate studies using data from general circulation models (GCMs).As examples, Hostetler and Clark (2000) used data from GENESIS (v.2.01) simulations and inferred ELA positions from mass balance profile curves for the Last Glacial Maximum on tropical glaciers.Rupper et al. (2009) also applied this kind of model to the region of Central Asia, using data from the NCEP/NCAR Reanalysis for the present and data from GCMs of phase 1 of PMIP for the Holocene.The aim was to reconcile Holocene glacier history with climate by quantifying the change in ELA for simulated changes in Holocene climate.
Details of the glacier mass balance model, which has been previously applied to Franz Josef Glacier, in New Zealand's Southern Alps, can be found in Anderson et al. (2006) and are briefly described here.This model calculates the mass balance gradient for any specific location, based on daily data of temperature and precipitation as a function of elevation.Elevation in the model is defined for each grid point from 0 C. Bravo et al.: Modelled glacier equilibrium line altitudes to 4000 m above sea level (a.s.l.) with steps of 20 m.For each elevation, the mass balance is calculated based on where ṁ is the mass balance rate, ċ the accumulation rate, and ȧ the ablation rate at time t and elevation z.
In glacier mass balance model, accumulation is defined as the portion of the daily precipitation that falls as snow when the daily average temperature is below a certain temperature threshold (T crit ).Previous studies have considered T crit being in the range of 0 to 2 • C (Radic and Hock, 2011).Therefore, water equivalent (w.e.) accumulation is calculated based on the daily information of mean temperature (T mean ) and total daily precipitation (p d ), calculated as In this case, T crit was assumed to be 1 • C (Anderson et al., 2006).
In the mid-latitudes, the ablation process is mainly controlled by melting and runoff (Rupper and Roe, 2008).Temperature is a good predictor of melt because incoming longwave radiation and turbulent heat fluxes are important terms in the energy balance that are closely related to air temperature (Ohmura, 2001;Oerlemans, 2001).The other major component of the energy balance, shortwave radiation, is also closely correlated with air temperature.
Ablation in the model is proportional to the mean daily temperature and occurs for values above 0 • C (Braithwaite, 1985;Hock, 2005).In this study, we calculated ablation using T mean when this is positive: where T + d is a positive daily temperature.Ablation is calculated by multiplying T + d by a factor that relates temperature and ablation, the degree day factor (DDF).The DDF (mm w.e.d −1 C −1 ) corresponds to the amount of melting (of ice and snow) per day, which occurs when temperatures are higher than 0 • C.This parameter shows great spatial variability and, in general, is higher for ice and lower for snow due to the high albedo of the latter that reduces the absorption of shortwave radiation (Braithwaite, 1995).In this study we use values of 6 and 3 mm w.e.d −1 C −1 for ice and snow, respectively.These values correspond to the same values used by Woo and Fitzharris (1992) for reconstructing the mass balance for Franz Josef Glacier in the Southern Alps.Therefore ablation is estimated by the following relationship: In this study, we use DDF snow when the snow depth is greater than zero and DDF ice when the snow depth is equal to zero.
Note that, in this study, we assume that temperatures below zero do not contribute to melting (Hock, 2003), and any potential contribution of sublimation to the total ablation is neglected because it is likely small compared to melting.
By applying this model at different elevations, we obtain a glacier mass balance curve (specific mass balance with altitude).The ELA occurs where the mass balance equals zero.For the purpose of this study we assumed that some parameters such as temperature and precipitation lapse rates, DDFs, and temperature threshold T crit are constant and equal for both the MH and the PI.Although this might not be strictly correct, our focus here is on the relative differences between the two periods rather than absolute values.
Given that the model has not previously been used in South America before, a preliminary validation was carried out by comparing model results with existing glacier mass balance data.Unfortunately, very few in situ glaciological mass balance measurements exist in southern South America (> 40 • S).The most recent published mass balance study in this area was at Mocho Glacier, located on the Mocho-Choshuenco volcano (39 • 55 S, 72 • 02 • W), which includes the hydrological year 2003/2004 (Rivera et al., 2005).
We assessed the performance of the model forced with climatological information coming from a regional climate model at 25 km resolution and forced by the ERA-40 reanalyses (PRECIS-ERA40) in Mocho Glacier in the Chilean Lake District (∼ 40 • S).Because the ERA-40 reanalyses cover the period 1957-2001, we cannot directly compare the same year analysed in Rivera et al. (2005) with the glaciological method (use of stakes for measurement of change in the ice/snow surface).We therefore run the mass balance model for the closest 10-year period (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) to the actual observed year in order to obtain an approximation to the climatological ELA.We use three precipitation lapse rates: 0, 0.00252, and 0.02 mm m −1 .This comparison (not shown) showed a reasonable agreement between the observed and modelled ELA, for the low precipitation lapse rates (0 and 0.00252 mm m −1 ), with a modelled ELA of about 100 to 150 m below the observed ELA, which is in agreement with the glacier recession documented for this glacier in the last 30 years (Rivera et al., 2005).

Model inputs: PMIP2
We use daily GCM temperature and precipitation values for computing the degree day model obtained from simulations carried out under PMIP2; see Braconnot et al. (2007) for model setup and boundary conditions.
Although PMIP is currently in its third phase (PMIP3), we used the modelling outputs of PMIP2 given that daily data were not available for the most recent phase when this study began.We analysed seven models of the PMIP2 initiative (Table 1).
We compared the PI outputs with gridded temperature and precipitation data from CRU TS3.10 and CRU TS3.10.01  (Harris et al., 2014), respectively, and with weather station observations to assess the climate model results (Figs.S1 to S4 in the Supplement).The glacier mass balance model was driven for 50 years in order to be able to capture the averaged influence of interannual variability, with daily temperature and precipitation data derived from the MH and PI experiments.With the purpose of comparing the ELA between the two periods, we calculated the mean mass balance for 50 years.Hence, the ELA calculated for each period corresponds to a long-term or climatic ELA.Temperature data were calculated for different elevations using a standard lapse rate of 6.5 • C km −1 .Due to the scarcity of available precipitation observations at high altitude, especially in Patagonia (Garreaud et al., 2013), precipitation was corrected using a regional mean observed gradient (in mm m −1 ) in both regions.The observed gradient was obtained using available latitudinal and altitudinal distributions of climate station data in both regions.Fitting the precipitation versus altitude distribution yielded a mean value of 0.00252 mm m −1 in Patagonia and 0.0038 mm m −1 in New Zealand.Given that the distribution of precipitation in mountainous regions is difficult to predict even under present-day conditions (Rowan et al., 2014), we use this simple approach to facilitate the process of mass balance modelling.In doing this, we are mindful that we are working at mountain range scale, and that the PMIP2 models do not represent the precipitation gradient very well, especially in the Southern Alps (Fig. S4).In addition, a constant precipitation factor (of 1.55) was also applied to account for the underestimation that lowresolution global models have for precipitation at high elevations (e.g.Rojas, 2006).
The results were averaged over six study zones.These zones correspond to the Chilean Lake District (CLD, approximately between 40 and 43 • S), Northern Patagonian Icefield (NPI, between 43 and 48 • S), Southern Patagonian Icefield (SPI, between 48 and 53 • S), and Cordillera Darwin (CD, between 54 and 55 • S) in South America, and the northern and southern sector of the South Island of New Zealand (NZN, between 41.5 and 43.5 • S, and NZS, 43.5 and 46 • S, respectively; Fig. 1).We also calculated climate differences between MH and PI over these six zones using monthly PMIP2 data and tested their significance using a t test in the case of temperature, and a non-parametric Wilcoxon-Mann-Whitney test in the case of precipitation with a significance level of 95 %.

Climate differences between the mid-Holocene and the pre-industrial period
Seasonal temperature differences between the MH and the PI (Fig. 2) are consistent between most models, with temperature anomalies of the same sign and small inter-model spread.Overall, in Patagonia, the models simulate colder conditions during the MH in austral summer (DJF), autumn (MAM), and winter (JJA), with average temperature anomalies of about −0.2 • C, about −0.5, and about −0.4 • C, respectively.During spring (SON) the PMIP2 models show temperatures 0.2 • C warmer.
In New Zealand (Fig. 2) the models show colder MH conditions in austral autumn (about −0.7 • C) and winter (about −0.4 • C), and warmer conditions in spring (about 0.3 • C).In summer the inter-model spread is larger, so that, on average, the temperature anomalies are not significant.In the annual mean, the temperature anomalies for South America and New Zealand are identical (about −0.2 • C).These temperature differences reflect the seasonal insolation difference between the two periods.Estimates of precipitation change show less consistency than for temperature (Fig. 3), and in several cases the models show precipitation anomalies of different sign within regions.Nevertheless, there are some regions and seasons for which the models show consistent precipitation changes.For example, during austral summer and autumn the models suggest that the climate was wetter during the MH compared to PI, in the CLD and the Patagonian icefields.In general all zones exhibit drier winters than the PI; spring was drier in the CLD and NPI, somewhat wetter in the SPI and CD, and marginally drier in the New Zealand zones.We find that the CLD was wetter in summer and autumn, experienced no change in winter, and was drier in spring.Note that none of the precipitation changes are statistically significant (p < 0.05, at 95% confidence interval).

ELA calculations and differences
For ELA calculations, we excluded the FOAM model due to its unsatisfactory simulation of the PI climate results (Figs.S1 to S4).The spatial distribution of the PI mean ELA based on six PMIP2 models in Patagonia (Fig. 4) shows that the ELA values are higher in the northern section of the study area, with maximum values above 2000 m a.s.l.(mean value of 1800 m a.s.l.) in the CLD.To the south, the ELA decreases gradually, reaching altitudes below 1000 m a.s.l. in SPI (mean value of 960 m a.s.l.) and CD (mean value of 840 m a.s.l.; see Table 2).Our results also show that, in general, the inter-model variability (1 standard deviation) of the ELA estimation is small.One exception is the inter-model variability observed in the SPI, where the maximum standard deviation is 250 m (mean value of 140 m), in a region with mean ELA of about 950 m a.s.l.
PI ELA estimates in New Zealand (Fig. 5) are higher in the northern part of the South Island (maximum values around 1800 m a.s.l.) and slowly decrease to values of 1400 m a.s.l. at the southern tip of the South Island (Table 2).These modelled values fall within the range of PI ELAs reconstructed from moraine evidence in the South Island (Lorrey et al., 2013).The inter-model spread evaluated with the standard deviation is small in the northern part (148 m).South of 43 • S, the inter-model spread becomes larger, with values between 150 and 180 m on the western flank and up to 200 m on the eastern flank of the Southern Alps.The mean value in this zone is 1530 m a.s.l.(Table 2).
As for the multi-model mean ELA differences, in Patagonia (Fig. 6a) and in the Southern Alps (Fig. 6b) the ELA was lower during the MH compared to PI; however the magnitude of change is relatively small: in Patagonia the mean difference is ∼ 20 m in all zones, and in the Southern Alps it is ∼ 30 m in both zones.Besides the small estimated ELA variations, it is important to highlight the consistency between ELA differences calculated by the PMIP2 models.In the Southern Alps, all of the six models indicate a negative sign in the ELA differences between the MH and the PI (Fig. 6).In Patagonia at least four models show negative differences between MH and the PI in almost the entire domain, with five models showing a lower ELA during the MH in some parts of the CLD and SPI zones and six models showing the same result in the west of the SPI zone (Fig. 6).

Difference between mid-Holocene and pre-industrial climates
Rojas and Moreno (2011) evaluated the full PMIP2 MH simulations for the climatic conditions in Patagonia and New Zealand.They found that both regions received less precipitation during a colder accumulation season, and more precipitation during a warmer ablation season.Therefore they suggested, on a qualitative basis, that the temperature and precipitation anomalies could effectively lead to larger glacier extents during the MH and hence explain neoglaciation in both regions.This paper goes a step further towards understanding the effects of climatic conditions on glaciers and neoglaciations, using those conditions to drive a glacier mass balance model.
The differences in temperatures between the MH and the PI found in this study are similar to those determined by Ackerley and Renwick (2010) for New Zealand, as well as Rojas and Moreno (2011) for South America and New Zealand.Both studies analysed data from PMIP2 models, but each used a different subset of models.Ackerley et al. ( 2013) use a regional simulation of higher spatial resolution to simulate the MH climate, and they also find a similar temperature pattern to that found in this study.Given that all these studies indicate that the MH was cooler than the PI in the autumn months and warmer in the spring months, we conclude that the temperature signals are robust across different subsets of PMIP2 simulations for the MH.
Our results indicate mostly wetter conditions during summer (DJF) and drier conditions in winter (JJA), in accordance with Rojas and Moreno (2011).For the autumn (MAM) and spring (SON) seasons there is dipole-like signal, with posi-  tive precipitation anomalies in the northern regions and drier conditions in the southern regions in MAM and the opposite for SON.These results are also in fair accordance with Rojas and Moreno (2011).In the New Zealand case, we find large inter-model spread between seasons, except during JJA, where we find clear drier conditions.Precipitation changes are slightly different than those shown in Ackerley and Renwick (2010), which in turn do not agree with the results of Rojas and Moreno (2011).In summary, we found small changes in precipitation and large inter-modal spread, so existing studies discussed here give slightly different results.

Differences between mid-Holocene and pre-industrial ELAs
We observe that the mass balance model applied to Patagonia and New Zealand is able to capture the expected differences in the climatological ELA associated with the climate conditions estimated for the MH and the PI.Our results show that, during the MH, the ELA could have been between 20 and 30 m lower than during the PI in Patagonia and New Zealand.
We propose that the results of the modelled ELA differences can be explained mainly by the significant and consistent differences in modelled temperatures observed.The impact of the precipitation anomalies is more difficult to assess, given that the climate data are heterogeneous.This suggestion is consistent with the idea that glaciers from midlatitudes are more sensitive to changes in temperature than to  changes in precipitation (Anderson and Mackintosh, 2006).Moreover, we suggest that the observed differences in climatological ELA are mainly driven by changes in the annual temperature cycle in these temperate regions.In Patagonia, ablation dominantly occurs between September and March (spring and summer months), whereas accumulation occurs from April to August (autumn and winter months; Rodbell et al., 2009).In New Zealand most of the ablation occurs between November and April (summer and parts of spring and autumn), whereas accumulation occurs from May to October (winter and parts of autumn and spring).However, both regions experience significant interannual variations, where accumulation or ablation sometimes persists for a longer period.
In Patagonia and New Zealand the lower summer temperatures observed during the MH imply less energy input and hence lower amounts of ablation.Although the opposite happens in spring, when the higher temperatures of the MH indicate greater ablation, we suggest that this change was not sufficient to balance the impact of the lower summer temperature on the mass balance.In addition, the lower temperatures observed during autumn and winter would increase the percentage of precipitation falling as snow rather than rain during the MH, as hypothetically suggested by Sagredo and Lowell (2012) and Rodbell et al. (2009).This is particularly critical in the Southern Alps, where a significant portion of the modern precipitation falls roughly at the elevation of the ELA.A temperature drop in this area would result an increase in the proportion of snow precipitation (Rother and Shulmeister, 2006).This can be particularly important in autumn and spring, when temperatures in the vicinity of the ELA are typically −1 to 2 • C. Additionally, precipitation during winter is higher during the MH in almost all the PMIP2 models in all zones (Fig. 3); this also contributes to accumulation and therefore a lower ELA in the MH with respect to PI.

ELA sensitivity to model parameters
We performed sensitivity runs to increase the robustness of the modelling results in the Patagonian and Southern Alps sectors.This motivated by the small differences in modelled ELAs and the lack of constraints on important parameters owing to the scarcity of measurements, especially in Patagonia.We investigated the sensitivity of ELA to the precipitation lapse rate and the DDF of snow and ice.
We assessed precipitation lapse-rate values of 0, 0.001, and 0.02 mm m −1 considering the observed precipitation gradient in Patagonia and the Southern Alps of New Zealand (0.00252 and 0.0038 mm m −1 , respectively).The sensitivity in the ELA difference to the precipitation lapse rate has a maximum of 15 m in the northern part of the South Island (Fig. 7; glaciers do not exist in this zone).At the latitude of the northern glaciers of the Southern Alps (approximately 43 • S) the sensitivity is close to 10 m.Sensitivity declines southward (see Fig. 7).In Patagonia (Fig. 8), the CLD has a maximum sensitivity of 6 m.This value is lower in the NPI zone (2 to 3 m) and close to 5 m in the SPI.From both figures it is clear that, in almost all of the study zone, for higher precipitation lapse-rate values, the ELA differences between the two periods become larger.We therefore conclude that the small ELA differences in Patagonia and the Southern Alps are significant and robust to this parameter and therefore the results presented are the most conservative modelled ELA differences.
We assessed DDF values of 4, 8, and 10 mm d −1 • C −1 for ice and 2 and 4 mm d −1 • C −1 for snow, within a range of the- oretical values used in the glacier modelling (3 mm d −1 • C −1 for snow and 6 mm d −1 • C −1 for ice).DDF sensitivity is even lower in the Southern Alps, with a maximum of 3 m in the northern part and also a reduction in the sensitivity to the south (Fig. 7).In Patagonia, sensitivity is 3 to 4 m along the Andes (Fig. 8).

Comparison of geomorphically reconstructed ELA and model results
In the following paragraphs we assess our estimates of ELA change against some records of neoglacial activity in both study areas In Lago General Carrera (named Lago Buenos Aires on the Argentinean side, 46 • 30 S), central Patagonia, it has been shown that glaciers advanced around 6200 years BP (Douglass et al., 2005).Geomorphic evidence at this site suggests that, during this glacial advance, the ELA dropped to 1100 m a.s.l., a 300 m difference with respect to the present position estimate for small isolated cirque glaciers.Further evidence of glacier activity is found by Harrison et al. (2012), who determined ages of 5700 years BP for a moraine located to the west of the NPI (46 • 36 S, 73 • 57 W) associated with the San Rafael Glacier.Recently, Strelin et al. (2014) found evidence for glacier advance between 6000 and 5000 years BP based on moraine dating on the east side of the SPI (Lago Argentino), specifically associated with the Upsala, Agassiz, and Frías glaciers.However, in the latter two studies, ELA differences were not calculated.In the Southern Alps, geomorphic and geochronological evidence suggest that the Tasman and Mueller glaciers (43 • 50 S, 170 • E) advanced around 6740 ± 160 years BP (Schaefer et al., 2009;age updated in Putnam et al., 2012).Putnam et al. (2012) demonstrated that a MH glacial advance also occurred at the Cameron Glacier (∼ 43 • 20 S, 171 • E) at 6890 ± 190 years BP, suggesting a regional event.In Putnam et al. (2012), ELA estimated for MH was ∼ 140 m lower than present, and 110 m lower than present 180 ± 48 years ago.This suggests a fairly modest change (∼ 30 m) between the MH and PI, consistent with the results of the ELA modelling in this work.While there is a systematic difference between the PI ELA calculated by the model and modern observations, it is clear that relative changes in ELA are very similar between our work and the estimates by Putnam et al. (2012).
The qualitative agreement in the direction of change between our modelling results and geomorphic studies in these regions, despite absolute differences that are significantly smaller (of the order of the tens of metres), leads us to conclude that the mass balance modelling accounts for some but not all of the climatic differences between this two periods.
Several caveats in our study account for more quantitative agreement in the absolute value of the ELA.(1) We compare the MH simulations with PI conditions, which are different from late 20th century climate.Reconstructions based on geologic evidence, on the other hand, are compared to late 20th century conditions.(2) Glacier advances during the mid-Holocene did not necessarily coincide with a precise time slice, namely 6000 years BP.(3) The spatial scale of individual glaciers and the coarse resolution of climate models hinder a direct comparison.(4) Important uncertainties are present in model parameters, especially those related to the spatial distribution of precipitation and DDFs.(5) The magnitude of the glacier expansion and mass balance from a given ELA change depends on local conditions and characteristics of the glaciers, for example glacier bed slope and hypsometry (Oerlemans, 2012;De Angelis, 2014).Glacier bed slope is a primary control on length sensitivity (Oerlemans, 2012); a glacier with a gentle bed slope, such as Upsala Glacier in the SPI, shows a high length sensitivity to ELA changes, estimated at ∼ −50 m of glacier length per metre of ELA increase (Oerlemans, 2012).With this in mind we expected large change in the accumulation and ablation areas even if the ELA oscillation is small (Mercer, 1965;Furbish and Andrews, 1984).In contrast, the steep Franz Josef Glacier shows a much smaller length sensitivity of ∼ −10 m of glacier length per metre of ELA increase.Glacier hypsometry is also an important control on the sensitivity of mass balance to change in ELA.Glaciers where the bulk of the area is located below the ELA are subject to the largest changes in mass balance for any given change in ELA (De Angelis, 2014).Considering all these aspects and limitations of the glacier mass balance model, we highlight agreement among both the sign of change and regional homogeneity within and between both study regions.This in-phase ELA response in Patagonia and New Zealand's South Island is also in agreement with glacier fluctuations observed during the 20th century, when glaciers seem to be in phase with similar climate forcing (Fitzharris et al., 2007).

Conclusions
A glacier mass balance model forced with PMIP2 simulations showed that southern mid-latitude glacier ELAs during the mid-Holocene (MH) were lower compared to preindustrial (PI) conditions.The robustness of these results is evaluated by using six different climate model data to run the mass balance model, as well as additional simulations varying a number of not well constrained parameters of the model such as the precipitation lapse rate and the snow and ice DDFs.The results of those sensitivity simulations showed that the ELA differences, although small, always had the same sign in New Zealand, i.e. lower ELAs during the MH compared to PI, and in most of the models in Patagonia.We are therefore confident that climatic conditions, as simulated by six PMIP2 models for MH conditions, could lead to larger glacier extents during this period.The main forcings of the modelled ELA differences are temperature differences.Significantly colder conditions during the summer, autumn, and winter months prevailed during the MH compared to the PI.In contrast, modelled precipitation changes were small and there was disagreement between models for the sign of change.Given that ELAs for the MH were consistently lower despite the range of precipitation data with which they were forced, our ELA results underline the evidence that temperate glaciers show a greater sensitivity to temperature changes (Anderson and Mackintosh, 2006).
Cooler temperatures during the MH would have affected glacier mass balance during both summer and winter.In summer and early autumn in the MH, less energy would be available for ablation, whereas from autumn to late winter, lower temperatures would cause a larger portion of precipitation to fall as snow, resulting in higher accumulation in the MH with respect to the PI, as well as a higher surface albedo, which reduces the amount of shortwave radiation available for melt.
This study provides new insights into understanding Southern Hemisphere mid-Holocene glacier conditions, demonstrating that orbital forcing inducing relatively coherent temperature changes is consistent with a hemispheric pattern of larger glacier extent at the MH compared to the PI period.
The Supplement related to this article is available online at doi:10.5194/cp-11-1575-2015-supplement.

Figure 1 .
Figure 1.Study area.(a) Schematic diagram showing some of the main climate characteristics in the southern mid-latitudes.(b) New Zealand South Island topography and the two zones of analysis (NZN: New Zealand north; NZS: New Zealand south).(c) Patagonia topography and the four zones of analysis (CLD: Chilean Lake District; NPI: Northern Patagonian Icefield; SPI: Southern Patagonian Icefield; CD: Cordillera Darwin).

Figure 4 .
Figure 4. Spatial distribution of pre-industrial equilibrium line altitude (ELA) in South America based on six PMIP2 simulations.(a) Mean ELA (m a.s.l.) and (b) inter-model variability of the ELA (1 standard deviation).White lines correspond to actual glacier extension according to the Randolph Glacier Inventory (RGI 3.2).

Figure 5 .
Figure 5. Spatial distribution of pre-industrial ELA in the South Island of New Zealand based on six PMIP2 simulations.(a) Mean ELA (m a.s.l.) and (b) inter-model variability of the ELA (one standard deviation).White lines correspond to actual glacier extension according to the Randolph Glacier Inventory (RGI 3.2).

Figure 6 .
Figure 6.Mid-Holocene minus pre-industrial ELA differences.Points indicate that the six models have a negative sign in the differences; asterisks indicate that five models have a negative sign and crosses indicate that four models have a negative sign.White lines indicate actual glacier outlines (Randolph Glacier Inventory 3.2).(a) Patagonia and (b) New Zealand.

Figure 7 .
Figure 7. Sensitivity of ELA differences in the Southern Alps to precipitation lapse rate and DDF of snow and ice.(a) Latitudinal topography profile of Patagonia (SRTM 1 km), (b) ELA difference sensitivity to precipitation lapse rate, and (c) ELA difference sensitivity to DDF of snow and ice.Grey and pink areas in (b) and (c) represent the range in ELA differences between MH and PI.

Figure 8 .
Figure 8. Sensitivity of ELA differences in the Patagonian Andes to precipitation lapse rate and DDF of snow and ice.(a) Latitudinal topography profile of Patagonia (SRTM 1 km), (b) ELA differences sensitivity to precipitation lapse rate, and (c) ELA difference sensitivity to DDF of snow and ice.Grey and pink areas in (b) and (c) represent the range in ELA differences between MH and PI.

Table 1 .
PMIP2 models used in this study.

Table 2 .
Mean values of ELA for each zone.