Inferring Past Climate from Moraine Evidence Using Glacier Modelling

Inferring past climate from moraine evidence using glacier modelling by Alice Marie Doughty A thesis submitted to Victoria University of Wellington in fulfilment of the requirements for the degree of Doctor of Philosophy in Physical Geography. Victoria University of Wellington 2013 ii Abstract Glacier length fluctuations reflect changes in climate, most notably temperature and precipitation. By this reasoning, moraines, which represent former glacier extent, can be used to estimate past climate. However, estimating palaeoclimate from moraines is not a straight-forward process and involves several assumptions. For example, recent studies have suggested that interannual stochastic variability in temperature in a steady-state climate can cause a glacier to experience kilometre-scale fluctuations. Such studies cast doubt on the usefulness of moraines as climate proxy indicators. Detailed glacial geomorphological maps and moraine chronologies have improved our understanding of the spatial and temporal extent of past glacial events in New Zealand. Palaeoclimate estimates associated with these moraines have thus-far come from simple methods, such as the accumulation area ratio, with unquantifiable uncertainties. I used a numerical modelling approach to approximate the present-day glacier mass balance pattern, which includes the effects of snow avalanching on glacier mass balance. I then used the models to reconstruct palaeoclimate for Lateglacial and Holocene glacial events in New Zealand, and to better understand moraine-glacier-climate relationships. The climate reconstructions come from simulating past glacier expansions to specific terminal moraines, but I also simulated glacier fluctuations in response to a previously derived temperature reconstruction, and to interannual stochastic variability in temperature. The purpose behind each simulation was to identify the drivers of significant glacier fluctuations. The modelling results support the hypothesis that New Zealand moraine records reflect past climate, especially changes in temperature. Lateglacial climate was reconstructed to be 2-3◦ C lower than the present day. This temperature range agrees well with previous estimates from moraines and other climate proxy records in New Zealand. Modelled temperature estimates for the Holocene moraines are slightly colder than those derived from simpler methods, due to a non-linear relationship found between snowline lowering and glacier length. This relationship results from the specific valley shape and glacier geometry, and is likely to occur in other, similarly-shaped glacier valleys. The simulations forced by interannual stochastic variability in temperature do not show significant (>300 m) fluctuations in the glacier terminus. Such fluctuations can not explain the Holocene moraine sequence that I examined, which extends >2 km beyond the present-day glacier terminus. Stochastic temperature change could, however, in part, cause fluctuations in glacier extent during an overall glacier recession. Modelling shows that it is also unlikely that glaciers advanced to Holocene and Lateglacial moraine positions as a result of precipitation changes alone. For these reasons, temperature changes are a necessary part of explaining past glacier extents, especially during the Lateglacial, and the moraines examined here likely reflect changes in mean climate in New Zealand. The glacier modelling studies indicate that simpler methods, such as the accumulation area ratio, can be used to appropriately reconstruct past climate from glacial evidence, as long as the glacier catchment has a straight forward geometry, shallow bed slope and no tributary glaciers. Non-linear relationships between climate change and glacier length develop when valley shape is more complex, and glaciers within these systems are probably better simulated using a modelling approach. Using a numerical modelling approach, it is also possible to gain a greater understanding of glacier response time, length sensitivities, and estimates of ice extent in valleys within the model domain where geomorphic evidence is not available. In this manner, numerical models can be used as a tool for iii understanding past climate and glacier sensitivity, thus improving the confidence in the palaeoclimate interpretations. iv Acknowledgments The primary source of funding came from the New Zealand Government under a New Zealand International Research Doctoral Scholarship. Additional funding was provided by Victoria University of Wellington (Antarctic Research Centre, Science Faculty, School of Geography, Environment and Earth Sciences) and the Comer Science and Education Foundation. R. Drews, T. Kerr, L. Kees, A. Mackintosh, and Mt. Hutt Helicopters assisted in field work on Cameron Glacier. Discussions with R. Dadic, H. Horgan, N. Golledge, H. Purdie, L. Carter, A. Lorrey, J. Renwick, R. Mucadam, and A. Tait were particularly helpful. I am grateful for the support provided by my supervisors, A. Mackintosh and B. Anderson. In addition to teaching me about glacier modelling, they encouraged me to tutor, attend conferences, seek additional funding, and organise fieldwork. Their advice has helped me grow from a student to a collaborator and I look forward to working with them on future research projects. I also appreciate the time I have spent with my supervisors’ families, who have helped me feel at home in a foreign country. Continued support and motivation comes from my colleagues from the University of Maine and Lamont-Doherty Earth Observatory. I want to thank G. Denton for reminding me how exciting the mysteries of the world are, A. Putnam for encouraging additional simulations to answer significant questions, D. Barrell for supplying beautiful maps and detailed comments, T. Chinn for enthusiastic discussions about New Zealand glaciers, and to the late B. Andersen, thank you for teaching me how to map Norwegian style. Mental sanity support was provided by M. Donaldson and my officemates, in particular M. Patterson, L. Kees, R. Rhodes, R. Gavey, and L. Cohen. My friends and flatmates (especially J. Latchford) as well as the staff at Ferg’s Kayaks were particularly helpful in motivating me to finish successfully. I thank my family for their unconditional support and love. I also appreciate the contributions, including artwork and story telling, from the newest member of my family, my niece Olivia (born 3 August, 2009). v vi Contents Abstract ii Acknowledgments v Table of Contents vi List of Figures xii List of Tables xix 1 2 Introduction 1 1.1 Former glaciers and climate change . . . . . . . . . . . . . . . . . . . . 1 1.2 Organisation of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Statement on the contributions made to this thesis by the author, supervisors, and collaborators . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Background 11 2.1 Moraines and palaeoclimate . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 New Zealand Lateglacial and Holocene palaeoclimate . . . . . . . . . . . 13 2.3 Lateglacial and Holocene glacial geology of New Zealand . . . . . . . . 15 2.3.1 Mapping Lateglacial and Holocene moraines . . . . . . . . . . . 16 2.3.2 Dating of moraines . . . . . . . . . . . . . . . . . . . . . . . . . 17 Climate of the Southern Alps . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4.1 Southern Alps climate variability . . . . . . . . . . . . . . . . . 19 2.4.2 Alpine meteorology . . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Glacier monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6 Glacier modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4 vii 3 4 Methodology 27 3.1 Modelling approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Energy-balance model . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Shortwave radiation . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 Albedo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.3 Longwave radiation . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2.4 Turbulent heat fluxes . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Gravitational snow mass transport and deposition model . . . . . . . . . 38 3.4 Models that describe glacier flow . . . . . . . . . . . . . . . . . . . . . . 43 3.4.1 Ice-flow model . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.5 Model applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Evaluation of Lateglacial temperatures in the Southern Alps of New Zealand based on glacier modelling at Irishman Stream, Ben Ohau Range 51 4.1 4.2 4.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.1.1 Glaciers as a climate proxy . . . . . . . . . . . . . . . . . . . . . 53 4.1.2 Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Modelling glacier extent . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.1 Input data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2.2 The energy-balance model . . . . . . . . . . . . . . . . . . . . . 58 4.2.3 The ice-flow model . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2.4 Steady-state simulations . . . . . . . . . . . . . . . . . . . . . . 63 4.2.5 Transient runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3.1 Experiment 1: Steady-state simulations of the ACR Irishman glacier extent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3.2 Experiment 2: Transient run driven by chironomid-derived temperature reconstructions . . . . . . . . . . . . . . . . . . . . . . viii 67 5 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Optimising simulated snow depth patterns on an avalanche-fed glacier in New Zealand using GPR and a gravitational mass transport model 75 5.1 5.2 5.3 5.4 6 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.1.1 Geographic Setting . . . . . . . . . . . . . . . . . . . . . . . . . 78 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2.1 Snow depth measurements . . . . . . . . . . . . . . . . . . . . . 80 5.2.2 Modelling snow distribution . . . . . . . . . . . . . . . . . . . . 82 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.3.1 GPR snow depth estimations . . . . . . . . . . . . . . . . . . . . 88 5.3.2 Ability of the model to replicate snow distribution . . . . . . . . 92 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 What do Holocene moraines at Cameron Glacier, New Zealand tell us about past climate change? 101 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.1.1 6.2 6.3 6.4 Site description . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.2.1 Bedrock and ice-surface topography . . . . . . . . . . . . . . . . 107 6.2.2 Climatic input data . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2.3 Mass balance model . . . . . . . . . . . . . . . . . . . . . . . . 109 6.2.4 2D ice-flow model . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.2.5 Experimental set-up . . . . . . . . . . . . . . . . . . . . . . . . 113 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.3.1 Steady-state experiments . . . . . . . . . . . . . . . . . . . . . . 118 6.3.2 Stochastic temperature forcing . . . . . . . . . . . . . . . . . . . 125 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 ix 6.5 7 6.4.1 Can stochastic variability in temperature cause the length of Cameron Glacier to fluctuate to its Holocene extents? . . . . . . . . . . . . 128 6.4.2 How do these results differ from previous studies of interannual variability influencing glacier length and why? . . . . . . . . . . 129 6.4.3 How do the modelled palaeoclimate estimates and ELA changes compare to previous estimates? . . . . . . . . . . . . . . . . . . 130 6.4.4 What do the Cameron Glacier Holocene moraines tell us about climate changes? . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.4.5 How do model results fit into the real world? . . . . . . . . . . . 133 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Synthesis 137 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.2 Contributions made in this thesis toward understanding palaeoclimate and glacier sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.3 Research questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.4 7.3.1 What magnitude of past temperature and precipitation change do Lateglacial- and Holocene-aged moraines in New Zealand represent? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 7.3.2 Are these estimates of past temperature and precipitation coherent with other climate proxy records in New Zealand? . . . . . . . . 140 7.3.3 To what extent does the re-mobilisation of snow by avalanches affect the present-day mass balance of New Zealand glaciers? . . 141 7.3.4 Does including a model of snow redistribution in a glacier model change the palaeoclimate reconstructions made? . . . . . . . . . 142 7.3.5 Could natural, stochastic variability in climate have caused Holocene glacier fluctuations? . . . . . . . . . . . . . . . . . . . . . . . . 143 7.3.6 What parameters are the modelled glaciers most sensitive to, and how does parameter choice affect the climate reconstructions? . . 143 7.3.7 How do the palaeoclimate estimates derived from modelling glaciers differ from previous estimates using simpler methods? . . . . . . 145 7.3.8 What can we learn from dynamic glacier models that we could not gain from using simpler methods? . . . . . . . . . . . . . . . 146 Proposed future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 x 7.4.1 Present-day climate and glaciology . . . . . . . . . . . . . . . . 148 7.4.2 Model experiments . . . . . . . . . . . . . . . . . . . . . . . . . 149 A Excerpt from: The anatomy of ‘long-term’ warming since 15 kyr ago in the Southern Alps of New Zealand based on net glacier snowline rise 153 A.1 Site description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 A.2 Excerpt from Kaplan et al., submitted . . . . . . . . . . . . . . . . . . . 155 A.2.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 A.2.2 Modeling and ELA Results . . . . . . . . . . . . . . . . . . . . . 155 A.3 Subsequent simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 B Ground penetrating radar profiles from the Cameron Glacier neve 161 B.1 GPR survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 B.1.1 Low-frequency data . . . . . . . . . . . . . . . . . . . . . . . . 161 B.1.2 High-frequency data . . . . . . . . . . . . . . . . . . . . . . . . 162 References 181 xi xii List of Figures 2.1 Map of shaded topography of central South Island, New Zealand with the last glacial maximum ice extent (pink outline Barrell et al. (2011)) and the Main Divide of the Southern Alps (dashed red line). Numbers mark locations of 1) Franz Josef Glacier, 2) Hokitika, 3) Arthurs Pass, 4) Christchurch, 5) Irishman Stream, 6) Whale Stream, and 7) Cameron valley. 16 2.2 Map of 30-year mean annual precipitation surface (Stuart, 2011), overlain on shaded topography of central South Island, New Zealand and the Main Divide of the Southern Alps (dashed black line). In the central Southern Alps, most of the annual precipitation falls west of the main divide due to a strong orographic effect. . . . . . . . . . . . . . . . . . . . . . . . . . 20 Modelled versus measured temperatures at the location of a climate station in Irishman basin. Modelled temperature values (daily interpolated temperature, T , smoothed with a thirty-day running mean) were interpolated from climate stations at <800 m asl over the 2010-2011 (a) and 2011-2012 (b) intervals. Measured temperature data (black line) are compared to modelled temperatures that were calculated using the 2010-2011 monthly varying lapse rate (red line), using a new variable lapse rate (purple line, panel b only) and a constant lapse rate (blue dashed). Using a constant lapse rate achieves a poor match, especially during winter months (Jun-Sep). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 One-dimensional example of snow redeposition from the MTD parameterisation modified from Gruber (2007). Panels show avalanche profiles resulting from changing parameter values, a) increasing deposition limit (Dlim ) leads to thicker, less-extensive deposits, b) increasing total incoming snow will increase the extent of the avalanche deposit, c) increasing slope limit (βlim ) causes snow to stay close to the valley wall, and d) increasing gamma (γ) also results in less-extensive, thicker deposits. For these tests, the total incoming snow is a 100 m s.w.e. column deposited at the circle labelled ‘input’, which is an unrealistic process but helps to highlight the shape of the deposit. . . . . . . . . . . . . . . . . . . . . . 39 Plots of slope angle versus maximum deposition to show parameter influences on the slope-deposition curve. Parameters evaluated are a) Dlim (yellow circles) set to 0.4 (dotted line), 0.6 (dashed line), and 0.8 (solid line) m s.w.e., b) γ = 0.075 (dotted line), 0.1 (dashed line), and 0.125 (solid line), and c) βlim (green circles) set to 35 (dotted line), 40 (dashed line), and 45◦ (solid line). . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Cameron Glacier model domain showing a) slope (◦ ) at 25 m grid resolution, b) slope at 50 m grid resolution, c) modelled mass balance (m a−1 ) using the MTD model, and d) modelled mass balance without the MTD model for the 2009-2010 simulation. The black line represents part of the end of summer snowline measured on 19 March, 2010 using a hand-held GPS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.1 3.2 3.3 3.4 xiii 4.1 4.2 4.3 4.4 4.5 Location of the study area in central South Island, New Zealand. The Ben Ohau Range is outlined in green and Irishman basin is in the white box. Boundary Stream tarn (BST) is represented by the white dot southeast of Irishman basin and the white triangles are the locations of Aoraki/Mt. Cook village (north of the Ben Ohau Range) and Twizel (southeast of Ben Ohau Range). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Modelled versus measured temperatures at the location of a climate station in Irishman basin. Modelled temperature values (daily interpolated temperature, T , smoothed with a thirty-day running mean) were interpolated from climate stations at <800 m asl over the 2010-2011 interval. Measurements (black line) are compared to modelled temperatures that were calculated using a monthly varying lapse rate (red line) and a constant lapse rate (blue dashed). Using a constant lapse rate achieves a poor match, especially during winter months (Jun-Sep). . . . . . . . . . . . . . 59 Modelled extent of Irishman glacier at 13 ka where ∆T -2.5◦ C and ∆P +11%. We show the (a) mapped geomorphology in the Irishman basin, for details see Kaplan et al. (2010), (b) modelled ice thickness overlaying the moraine map, (c) modelled mass balance with modelled ELA (green dashed line), and (d) previous estimates of Irishman glacier extent (K1 blue, K2 red) and ELA (K1ELA blue dashed, K2ELA red dashed) (Kaplan et al., 2010) based on geometric and AAR reconstructions, compared to this study (green). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Experiment 1 results - Inferred ∆P -∆T reconstructions showing the sensitivity of the climate reconstruction to chosen parameter values. Curves identical to the ‘Optimal’ curve are not shown (I, Uc , A, and Zice =0.002 m). The x-intercepts of the different curves are reported in Table 4.3. . . . . . 70 Experiment 2 results - BST PLS-smooth temperature reconstructions from Vandergoes et al. (2008) (black) from 15 to 11 ka with modelled glacier length (blue). Orange and yellow diamonds represent the moraine positions dated to 13.0±0.5 and 12.1±0.5 ka respectively. Glacier length was calculated by identifying the ice with a thickness >10 m farthest from the headwall. Plots in the right panel show a) PLS-raw, b) WAPLS-raw where tick marks show the sampling resolution, and c) WAPLS-smooth inferred summer mean ∆T from 15 to 11 ka. Glacier lengths resulting from these three reconstructions extended beyond the 13 ka moraine and are not shown here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 xiv 5.1 5.2 5.3 Topographic map of the Cameron Glacier catchment. a) An oblique aerial photograph of the Cameron Glacier showing the steep valley walls, partially debris-covered terminus, and end of summer snowline (red dotted line). The black line delineates the ridgeline of the upper Cameron catchment. Photograph taken 5 March, 2010 by A. Willsman. b) Slope values for a section of the model domain (25 m grid resolution) with labelled GPR profiles on Cameron Glacier. Note the large area on the northern slopes that are >60◦ (black cells). Each profile is ∼200 m long. Profile colours are for differentiating lines only. c) Topographic map with elevation contours, ridgelines (dotted-dashed red line), glacierised area, and mountain peaks (Mt. Arrowsmith (2781 m asl) toward the southwest, and Jagged Peak (2705 m asl) to the northeast). The black star marks the location of the snow pit and the green star marks a location that does not receive avalanche accumulation. d) Shaded relief map of the central Southern Alps showing the location of study area (yellow square) and the location of the Mt. Potts weather station (white dot). Prominent glaciers, lakes, and rivers are also shown. . . . . . . . . . . . . . . . . . . . . . . 79 Slope angle (β) versus maximum snow deposition (Dmax ) showing the influence of parameter values on the shape of the Dmax -β curve. Parameters evaluated are a) Dlim = 0.4, 0.6, and 0.8 m, b) γ = 0.075, 0.1, and 0.125, and c) βlim = 35, 40, and 45◦ . Yellow circles mark Dlim values, green circles mark βlim values and the solid, dashed, and dotted curves relate to the different values used for each parameter in each plot. . . . . . 83 Interpreted snow depth from the December, 2009 Cameron Glacier lowfrequency GPR profiles. Profiles F-G, and Q-Q’ show cross sections of the glacier while profile E-L follows the central flowline (see key and Figure 5.1b). Note the strong reflector, interpreted as the end of summer glacier surface for 2009 (red dashed line) and deeper reflectors representing previous summer surfaces (not marked). Reflections in the E-L profile show old accumulated snow deposits below the 2009 summer surface. . . 89 5.4 Snow depth over the model domain (a-c) and along GPR profiles (d-f). Snow depth (m s.w.e.) from the a) optimal parameter set, b) less appropriate parameter set, c) model run without the MTD model, d) GPR data, e) optimal parameter set, and f) less appropriate parameter set. The snow distribution in the optimal model run reflects the observed accumulation pattern at Cameron Glacier better than the other runs. Panel c) also shows GPR lines for reference, a black star (white outline) to mark the location of the snow pit and a green star to mark a point of interest in the discussion. 90 5.5 Matrices showing relative comparisons (by colour) between modelled data and measured snow depth (a) gradient, (b) mean, and (c) range for each profile for each run. Red-coloured cells show where the model over estimated values, blue represents under-estimated values, and white is within the accepted range where modelled and measured are similar. Black arrows point to the ‘best-fit’ simulation and the final column (test 28) is the simulation without the MTD model. Profile key to the left is for reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 91 5.6 Graph showing the number of ‘acceptable fits’ (white cells in Figure 5.5 and Table 5.1) for the gradient (◦), mean (+), and range (∗) relative comparisons for each model run, with the total sum of the white boxes in bold. Note the overall decrease in the total (bold line) sum of white cells from MR1 - MR28. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.1 a) Oblique aerial photograph of the Cameron valley Holocene moraine sequence. This photograph shows the Cameron Glacier snout, the long, sharp-crested lateral moraines, and three of the targeted Holocene terminal moraines (M1 - M3). The most extensive moraine visible in this image (M3), located 3 km from the terminus, is dated to 8,190±230 years before 1950. Scree slopes are abundant in this catchment, as seen on the valley flanks, and the Cameron River has eroded through parts of the terminal moraines. View is towards the northwest, photo by T. Chinn. b) Oblique aerial photograph of Cameron Glacier with a dotted red line marking the end of summer snowline and black line marking the upper Cameron catchment ridgeline. This photograph also shows the upper Cameron catchment with its steep valley walls, which are prone to snow and rock avalanching. Photograph taken 3 March, 2010 by A. Willsman, view is towards the northeast. . . . . . . . . . . . . . . . . . . . . . . . . 103 6.2 Glacial geomorphological map of the southeast slopes of the Arrowsmith Range, showing the moraine sequences of the Cameron (northeast) and Ashburton (southwest) valleys (Barrell et al., 2011). The black line represents the central flowline from which glacier lengths were measured, and M1 - M5 mark each of the targeted moraine positions. Note the similarities in moraine sequences between the two valleys, although the Cameron sequence contains a larger number of ridges. Location diagram (top left) shows shaded topography (based on a DEM) of the central Southern Alps including the main drainage divide (red dashed line), prominent lakes, rivers, and glacierised area. The Arrowsmith Range is marked by a white box. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.3 White noise temperature forcing series with a standard deviation (σT ) of 1◦ C (blue diamonds) and the 30-year running mean temperature (red line). Note the relatively high peak in running mean around model year 1150, showing that the random temperatures can generate a small persistence in climate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.4 Glacier ice thickness and mass balance for the five different ice extents with no precipitation change. Note the simulated connection of tributary glaciers to the main trunk glacier at different extents of Cameron Glacier. Yellow arrows point to areas where the modelled ice extent and moraine map differ (see text). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.5 Glacier outline and ELA for each moraine extent. Panel a) glacier area extents overlapped for direct comparisons, and b) glacier outline and corresponding ELA (Table 6.3). . . . . . . . . . . . . . . . . . . . . . . . . 120 xvi 6.6 Temperature and precipitation change combinations needed to reach each of the moraines under 12 different parameter combinations (Table 6.2). Note that the palaeoclimate estimates overlap, at least slightly, with those produced for the neighbouring moraine. . . . . . . . . . . . . . . . . . . 122 6.7 Glacier length profile (m) response to a ±1◦ C temperature change from the steady-state temperature for M1 - M4. Note the greater length change caused by a -1◦ C temperature change, especially for M3 and M4. This asymmetric behaviour results from glacier geometry and valley shape. . . 123 6.8 Glacier length from M1 (equivalent to 4.3 km glacier length) plotted against ELA lowering from modern 2120 m asl from all tested steadystate model runs. Results from Putnam et al. (2012) based on AAR reconstructions are shown for comparison. Notice that the relationship between glacier length change and ELA change is linear for the late Holocene glacier extents, but differs for early Holocene extents, suggesting a strong topographic influence in this valley. . . . . . . . . . . . . . . . . . . . . 124 6.9 Glacier response to a white-noise forcing in temperature. Panel a) interannual temperature change for each of the transient experiments. Note the different amplitude in temperature change for the different tests (caused by σT ). Panel b) total ice volume in the model domain over time. Panel c) simulated Cameron Glacier length change for each of the stochastic variability experiments. Note the greater amplitude in glacier length fluctuations with increased σT (WNF-D, σT =2◦ C). . . . . . . . . . . . . . . 127 A.1 Shaded topography of the Ben Ohau Range (outlined in green) showing the locations of Irishman Stream and Whale Stream (white boxes). . . . . 154 A.2 Modeled ice extent for Whale Stream Glacier. (A) Various glacier outlines and ELAs using the AAR method, (B) modeled glacier extent showing ice surface elevation (20 m contours), accumulation (blue) and ablation (red) areas and ELA overlying the AAR method glacier extent, (C) modeled mass balance gradient ranging from -2 to +1 m a−1 , and (D) temperature-precipitation change curves with uncertainty for the Whale Stream Late Glacial ice extent. . . . . . . . . . . . . . . . . . . . . . . . 157 A.3 Parameter testing curves from the temperature-precipitation change combination needed for a Late Glacial Whale Stream ice extent. . . . . . . . . 158 A.4 Whale glacier simulated at 25 m resolution. The panels show a) ice thickness (m), b) mass balance (m a−1 ), c) slope at 25 m resolution (◦ ). . . . . 160 xvii B.1 Map of the Cameron Glacier. a) Topographic map with elevation contours, ridgelines (dotted-dashed red line), glacierised area, and mountain peaks (Mt. Arrowsmith (2781 m asl) toward the southwest, and Jagged Peak (2705 m asl) to the northeast). b) Shaded topography based on the DEM showing the location of study area (yellow square) in South Island, New Zealand. Prominent glaciers, lakes, and rivers are also shown. c) Slope values for a section of the model domain (25 m grid resolution) with labelled GPR profiles on Cameron Glacier. Profile colours are for differentiating lines only. . . . . . . . . . . . . . . . . . . . . . . . . . . 162 B.2 Plot of the GPS data for each profile. Each profile is labelled with its associated number, an arrow indicating the direction of travel from start to finish, and a colour to differentiate it from the other profiles. . . . . . . 164 B.3 Profiles 241, 243, and 244. . . . . . . . . . . . . . . . . . . . . . . . . . 165 B.4 Profile 245. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 B.5 Profiles 246, 247, 248, and 249. . . . . . . . . . . . . . . . . . . . . . . 167 B.6 Profile 250. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 B.7 Profile 251. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 B.8 Profile 252. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 B.9 Profile 253. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 B.10 Profile 254. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 B.11 Profile 255. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 B.12 Profile 256. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 B.13 Profile 257. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 B.14 Profile 258. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 B.15 Profile 259. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 B.16 Profile 260. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 B.17 Profile 261. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 B.18 Profile 262. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 xviii List of Tables 2.1 Climate information from Hokitika, Arthurs Pass, and Christchurch stations. Statistics from daily mean precipitation and temperature data from 1980 to 2011 (NIWA, Retrieved 2009-2011). Hokitika (west coast) and Christchurch (east coast) data are presented for near sea-level west to east comparisons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Energy balance model parameters, values, and references. . . . . . . . . . 38 3.2 Mass transport and deposition model parameter setup. . . . . . . . . . . . 40 3.3 Ice-flow model parameter options and the range of published values for Glen’s flow law coefficient. . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4 Model information summary for comparison between the different studies. 48 4.1 Measured monthly mean temperatures (Tmeas ) and calculated temperature ) for each month based on Irishman basin climate station lapse rates ( dT dz data (2010-2011, 43◦ 59’30”S, 170◦ 02’57”E, 2010 m asl). Annual mean measured temperature is ±1.2◦ C (assuming a value of 6.2◦ C for February) and the annual mean lapse rate is −5.4◦ C km−1 . . . . . . . . . . . . 59 Model parameter values. These values define our optimal parameter set, meaning they come from previously published studies optimal for modelling and/or are specific to the site location. Our sensitivity tests use these values, changing one value at a time to understand each parameter’s influence on our climate reconstruction. . . . . . . . . . . . . . . . . . . 63 Parameter name, symbol, value, and the ∆T needed to grow the glacier to the 13 ka moraine without a change in precipitation (x-intercept), with the optimal parameter set yielding a ∆T of −2.7◦ C. . . . . . . . . . . . . 67 Model run number with parameter combination and the number of profiles with an acceptable fit (where the relative comparison is between 32 and 23 ). The number of ‘acceptable fits’ per test for the gradient, mean, and range relative comparisons are listed as well as the total sum of ‘acceptable fits’ for each model run. MR28 is without the avalanche model. Note that model runs with βlim =35◦ perform relatively well (bold font). . 93 4.2 4.3 5.1 6.1 Cameron Glacier Holocene moraine information modified from Putnam et al. (2012). Information includes cosmogenic exposure ages (M1 - M4), associated age for M5 (Putnam et al., 2010a), ∆ ELA is relative to the AD1995 ELA for Douglas Glacier (2120 m asl). . . . . . . . . . . . . . . 105 xix 6.2 Optimal and tested model parameter values. The values with references come from previously published studies, whereas those with ‘S#’ are values tested in the sensitivity simulations. Names in bold text are the parameters we varied to define our palaeoclimate uncertainty. Our sensitivity tests involve changing one value at a time to quantify the influence each parameter has on our palaeoclimate reconstruction. . . . . . . . . . 114 6.3 Palaeoclimate estimates for the five targeted moraine positions. Temperature change range reflects ∆P values. Cameron Glacier length, maximum ablation, maximum accumulation, ELA, area, and AAR are offered for each moraine simulation. AAR values for the Ashburton are also offered, based on the modelled mass balance. . . . . . . . . . . . . . . . . . . . . 120 6.4 Sensitivity of modelled glacier to interannual temperature variability. Starting from the M1 glacier extent (WNF-A to WNF-D) where baseline/equilibrium glacier length is 4200 m and from the early Holocene moraine (WNF-E) where baseline glacier length is 5700 m. Below are values of modelled glacier length, including the percentage of time the glacier length spent behind or beyond its equilibrium length and the calculated standard deviation and skewness of the modelled glacier length variations. 126 6.5 Values describing the linear relationship between a change in ELA and its corresponding change in glacier length (y = ax + b). ‘Tributary glaciers’ refer to the Douglas and Marquee tributary glaciers specifically. . . . . . 130 7.1 Palaeoclimate estimates from different aged moraines for the three targeted locations. Ages are in 1000 years before present (ka) and the range of temperatures represents a precipitation change of -20%, 0%, and +20% relative to the present-day value. . . . . . . . . . . . . . . . . . . . . . . 139 A.1 Parameter names, values, and the temperature change needed to grow the glacier without a change in precipitation (x-intercept), with the ‘best fit’ model parameters yielding -2.1◦ C. . . . . . . . . . . . . . . . . . . . . . 156 A.2 Palaeoclimate estimates using two different modelling approaches. . . . . 159 B.1 Profile information. Direction codes are in reference to the true right or true left sides of the glacier or elevation, left to right (L2R), right to left (R2L), up to down glacier (U2D), down to up glacier (D2U). Profiles used to test the avalanche model are correlated with their alias. . . . . . . . . . 163 xx Chapter 1 Introduction 1.1 Former glaciers and climate change Glacier length fluctuations reflect changes in climate, most notably temperature and precipitation (Oerlemans, 2001, 2005). It follows that moraines, which are deposited by a glacier at its margins, can be used to reconstruct past glacier extent and therefore past climate. Previous studies have shown that palaeo-glacier extents can be simulated using numerical models. Models also afford the possibility of testing different hypotheses regarding past climate and/or glacier behaviour (e.g. Plummer and Phillips, 2003; Kessler et al., 2006; Laabs et al., 2006). Here, I use a modelling approach to help understand what moraines can tell us about palaeoclimate in New Zealand. In the process, I also investigate glacier mass balance processes, response times, and glacier length-valley shape relationships. Southern Hemisphere moraine sequences offer great potential for helping us to understand the regional and global climate system. The climate in the southern middle-latitudes is ocean-dominated, with only the Andes of South America and Southern Alps of New Zealand extending deeply into the region of the westerly winds. These mountain chains support present-day glaciers and show evidence for past glaciation (Gair, 1967; Porter, 1975a; Suggate, 1990; Hubbard et al., 2005; Ackert et al., 2008; Barrell, 2011; Kaplan et al., 2011; Strelin et al., 2011). The climate of the South Island of New Zealand is influenced by both sub-tropical and sub-polar, atmospheric and oceanic circulation. If a circulation type persists, it can lead to a shift in the climate of the South Island, influencing glacier mass balance (Lamont et al., 1999; Clare et al., 2002; Chinn et al., 2012). Comparing the timing of glacier advances, as determined by moraine ages, between the northern and southern middle latitudes offers clues as to the possible mechanisms of past climate change (Denton and Hendy, 1994; Schaefer et al., 2009; Kaplan et al., 2010; Putnam et al., 2010a). 1 The recently published central South Island glacial geomorphology map (Barrell et al., 2011) and emerging cosmogenic exposure age chronologies (Kaplan et al., 2010; Putnam et al., 2010a; Kaplan et al., submitted; Putnam et al., 2012) provide the focus and incentive to further examine moraine-glacier-climate relationships. The chronologies in these publications focus on Lateglacial (15,000 to 11,500 years ago) and Holocene (last 11,500 years) glacier fluctuations in several locations in the Southern Alps. Moraine ages are thought to represent the end point of a period in which a glacier terminus remained stationary (Schaefer et al., 2009). The significance of a moraine sequence, in relation to climate, should be better understood, and this requires a closer look at glaciological processes and glacier responses to climate within individual basins. In this study, I investigate Holocene and Lateglacial climates by modelling glacier extents in three different valleys. These valleys contain published moraine chronologies and detailed geomorphological maps, which allowed me to target dated moraine positions to estimate past climates. The majority of previous palaeoclimate estimates from glacial records in New Zealand are based on an accumulation-area ratio (AAR) approach (Porter, 1975b; Kaplan et al., 2010; Putnam et al., 2012), which assumes a simple and direct relationship between glacier size and equilibrium-line altitude (ELA). The change in ELA is then used to estimate past temperature, assuming a temperature lapse rate. I compare the published AAR palaeoclimate estimates to my simulated glaciers to test whether this simpler approach is appropriate. I also address several questions regarding glacier-climate relationships that could not have been answered with an AAR approach. In Chapter 4, the modelled glacier is forced by a temperature reconstruction from a palaeo-ecological proxy. In this type of study, the model is accounting for the response time of the glacier and how it reacts to the magnitude and duration of a temperature change in the climate proxy record. In Chapter 6, I examined the potential of a glacier to advance to its Holocene moraine positions when forced by natural year-to-year stochastic variability in climate rather than a systematic shift in climate. These transient model runs provide insight into glacier length reactions, which could potentially result in moraine deposition. 2 A numerical modelling approach, informed by present-day glaciology, can be used to test assumptions that are often made whilst interpreting glacier length changes. The model used here simulates glacier extent by estimating surface energy balance and iceflow within a model domain. Energy balance components are solved for explicitly, and thus topographic influences (e.g. aspect, shading, and elevation) are included. The modelled glacier extent is thus a product of the choice of parameter values as much as the prescribed climate change. Varying parameter values shows the sensitivity of the palaeoclimate reconstruction to each parameter and provides an indication of uncertainty for a reconstruction. Estimating palaeoclimate from moraines is not a straight-forward process and has several assumptions. The following principles form the basis of this study; (a) in steady-state tests, I assume glaciers deposited moraines while being in equilibrium with the climate, (b) in all model tests, glacier length changes are primarily related to changes in temperature and precipitation, (c) many of the present-day processes and climate variables (e.g. relative humidity, cloudiness, wind speed, and seasonality) are assumed to be the same for palaeo-simulations, and (d) lake calving and extensive debris cover has not influenced the glaciers that I have focused on to any large degree. In transient simulations, assumption (a) does not apply because we are interested in glacier length fluctuations, for example, in Chapters 4 (where modelled climate and the resulting glaciers are driven by a nearby palaeoclimate proxy) and 6 (where interannual temperature variability within a steady climate is considered). Assumption (b) is examined in the sensitivity tests of Chapters 4 and 6 where changing model parameter values could account for some of the glacier length change. Assumption (c) was not tested here because reliable climate proxy records do not yet exist for past relative humidity, cloudiness, or wind speed (Knudson et al., 2011). We assume that these climate variables are comparable to their present-day values. If these variables changed in the past, as might be expected, their influence on the modelled ice extent is less than temperature and precipitation changes, meaning that they are unlikely to greatly impact our results. The error bars associated with the temperature and precipitation change estimates somewhat account for possible errors introduced by making these assumptions. Lastly, assumption (d) is not examined 3 because geomorphic evidence does not suggest there were large lakes or rock avalanches in the valleys chosen. It is possible that simulating past glacier extents would require less of a cooling or precipitation increase if debris cover on the ablation area were included. I do not examine the effects of debris cover on past glacier extents because debris cover is unknown for the past and its effects are beyond the scope of this study. I aim to understand glacier-climate interactions using the following approach: 1. In the first experiment, I use a coupled energy-balance model (EBM) and ice-flow model to reconstruct past climate in the Irishman basin from a small cirque glacier in the central Southern Alps where a well-dated Lateglacial (15,000-11,500 years ago) moraine sequence exists (Kaplan et al., 2010). First, a suite of steady-state model tests provided an envelope of temperature and precipitation change combinations that resulted in ice extending to a 13,000 year old moraine. These results, along with model parameter sensitivity tests, show the relative importance of temperature and precipitation change during the Lateglacial. Second, timedependent simulations forced by a chironomid-derived temperature reconstruction from a nearby tarn are used to compare the resulting glacier length changes to the moraine age and position. This multiple-climate-proxy approach provides an interesting test of the coherence of palaeoclimate estimates in the central Southern Alps. 2. In the second experiment, I use a gravitational snow mass transport and deposition model (Gruber, 2007) to simulate measured snow depth at Cameron Glacier in the central Southern Alps. Different parameter value combinations were then used to optimise a parameter set that provided the best match between measured and modelled snow distribution. This model, which is part of the mass balance model used here, improves predictions of present-day snow deposition patterns on glaciers that receive snow avalanche accumulation (Machguth et al., 2006). Avalanche processes appear to be a significant source for accumulation on the glacier studied in this experiment, as well as many other glaciers in the Southern Alps, and their mass balance processes should be better understood. Modelled mass balance that includes 4 the avalanche model results in relatively low AAR values, which makes for an interesting comparison to more typical AAR values. These accumulation processes were likely to have occurred in the past as well. Thus, in an attempt to properly simulate past glacier mass balance, I used the avalanche model in simulations of past glacier extent (Chapter 6). 3. In the third experiment, I use an ice-flow model, coupled to an EBM, which includes an avalanche model, to simulate glacier fluctuations driven by steady-state changes in climate (centennial-scale), as well as stochastic temperature variability (interannual). Steady-state temperature and precipitation change combinations are calculated for four Holocene moraines and one Lateglacial moraine in the Cameron valley in the central Southern Alps. Several modelling studies have shown kilometrescale glacier fluctuations are possible in a steady climate in response to interannual temperature variability alone. This study is designed to better understand glacierclimate links and whether glacier fluctuations caused by stochastic temperature variability in a steady climate could cause the glacier to advance to its Holocene moraine positions. 4. In the final experiment, I reconstruct past climate for Lateglacial moraines in the Whale Stream catchment in the Southern Alps. This particular experiment is part of a larger study including glacial geomorphological maps and cosmogenic nuclidederived exposure-ages from moraines, neither of which formed a part of this thesis. The simulated glacier extent, ice thickness, mass balance, and associated AAR offer insight into the uncertainty associated with using a generic AAR (e.g. 0.65) in a palaeoclimate reconstruction. The simulations also complement the palaeoclimate estimates from the other modelled Lateglacial ice extents in this thesis (Chapters 4 and 6). Additional simulations include a finer horizontal grid spacing, an avalanche model, a higher characteristic sliding velocity, and daily temperature variability to show how differences in model parameter values affect the palaeoclimate interpretations. Collectively, by carrying out these four studies I have made the following contributions: 5 • New estimates of past climate for the Lateglacial and Holocene in New Zealand. • Comparisons between the glacial geology and modelled glacier fluctuations forced by a chironomid-based climate reconstruction. • Improved accumulation estimates for an avalanche-fed mountain glacier by including an avalanche parameterisation. • Comparisons between glacial geology and modelled glacier fluctuations forced by interannual stochastic variability in temperature. • Evaluation of the relationship between snowline lowering and glacier length change in areas with height-mass balance and hypsometric feedbacks. In all of these modelling studies, a common finding is that climate change, in particular temperature change, is the primary driver of glacier terminus fluctuations and subsequent moraine records. 1.2 Organisation of the Thesis The main chapters (Chapters 4-6) of this thesis are written as manuscripts that have either been published, or will be submitted to international journals. The background (Chapter 2), methodology (Chapter 3), and synthesis (Chapter 7) chapters, as well as the reference section, are formatted in traditional thesis style. As a consequence of this structure, repetition does occur, particularly in the description of the study sites and models. Chapters 4-6 each address one of the primary research objectives outlined in Section 2.7. Appendix A is unlike the main chapters in that it contains only an excerpt from a manuscript lead by M. Kaplan and others. The manuscript contains glacial geomorphological maps and cosmogenic exposure ages from moraines, neither of which were a result of this thesis and are therefore omitted here. Appendix B contains ground-penetrating radar (GPR) data used in Chapters 5 and 6. These separate chapters together represent a coherent attempt to understand the climatic significance of moraine records in the Southern Alps of New Zealand. 6 1.3 Statement on the contributions made to this thesis by the author, supervisors, and collaborators The design of this study, data acquisition, its interpretation and presentation, and writeup represent the efforts of the author. The ice-flow and energy-balance models were originally coded in MATLAB by B. Anderson (Victoria University of Wellington) and the gravitational snow mass transport and distribution model was translated into MATLAB from IDL by myself with help from R. Dadic and H. Horgan. Ground penetrating radar surveys were collected with the help of L. Kees (Victoria University of Wellington), T. Kerr (National Institute of Water and Atmospheric Research (NIWA)), and A. Mackintosh (Victoria University of Wellington). A. Mackintosh and B. Anderson, my supervisors, appear as co-authors on all the manuscripts. While their contributions to this thesis have been significant, their input to the papers that appear in the thesis was no greater than would be the case for a Ph.D. thesis following traditional protocol. The contributions of B. Anderson and A. Mackintosh were threefold: (i) provision of research funds; (ii) discussion and contribution of scientific ideas and concepts; and (iii) editorial input. Below is a list of manuscripts as they appear in the thesis. The list includes the co-authors, the journal for publication, and a brief account of the contributions of the various co-authors. Chapter 4: Doughty, A.M., Anderson, B.M., Mackintosh, A.N., Kaplan, M.R., Vandergoes, M.J., Barrell, D.J.A., Denton, G.H., Schaefer, J.M., Chinn, T.J.H., Putnam, A.E. (in press) Evaluation of Lateglacial temperature in the Southern Alps of New Zealand based on glacier modelling at Irishman Stream, Ben Ohau Range. Quaternary Science Reviews. DOI: 10.1016/j.quascirev.2012.09.013. I carried out the modelling simulations. Kaplan, Barrell, Denton, Schaefer, Chinn, Putnam, and I provided published cosmogenic ages and a detailed geomorphology map of Irishman basin. Vandergoes supplied published temperature reconstructions derived from 7 chironomid assemblages from Boundary Stream tarn. Chapter 5: Doughty, A.M., Dadic, R., Anderson, B.M., Mackintosh, A.N., Kerr, T., Kees, L., (in prep) Optimising simulated snow depth patterns on an avalanche-fed glacier in New Zealand using GPR and a gravitational mass transport model. Dadic and Anderson provided source code, support, and aided in experiment development. Kerr, Kees, and Mackintosh assisted in the fieldwork. Kees and H. Horgan (not listed) assisted with initial ground penetrating radar post processing. Chapter 6: Doughty, A.M., Mackintosh, A.N., Anderson, B.M., Putnam, A.E., Barrell, D.J.A., Denton, G.H., Schaefer, J.M. (in prep) What can the Holocene moraines at Cameron Glacier, New Zealand tell us about past climate change? I carried out the model simulations. Putnam, Barrell, Denton, Schaefer, and I provided cosmogenic ages and a detailed geomorphologic map of Cameron valley. Appendix A, excerpt from: Kaplan, M.R., Schaefer, J.M., Denton, G.H., Doughty, A.M., Barrell, D.J.A., Chinn, T.J.H., Putnam, A.E., Andersen, B.G., Mackintosh, A., Finkel, R.C., Schwartz, R. and Anderson, B. (submitted to Geology) The anatomy of ‘long-term’ warming since 15 kyr ago in the Southern Alps of New Zealand based on net glacier snowline rise. Here I provide an introduction to the setting described in Kaplan et al. (submitted) and present only the modelling description and results that were included in a publication lead by Kaplan. The publication focuses on the moraine mapping and dating efforts. I added a modelling component to provide estimates of palaeoclimate and mass balance, which allowed me to calculate an AAR and ELA for the targeted glacier extent. The ELA calculated from the model was then compared to the estimated ELAs, which was calculated from the accumulation area ratio method. 8 Appendix B GPR profiles used in Chapters 5 and 6 to calculate snow depth and ice thickness. Kerr, Kees, and Mackintosh assisted in the fieldwork. Kees and Horgan assisted with ground penetrating radar post processing. The raw data from this survey are available on the Victoria University of Wellington ‘tuawe’ network drive. 9 10 Chapter 2 Background Numerical modelling provides a potentially useful, physically-based method to examine glacier-climate interactions. In recent years, a number of developments have made modelling of past glacier extent in the Southern Alps of New Zealand more feasible. First, recently published glacial geomorphology maps and surface exposure dating (SED) ages of moraines in the central Southern Alps provide the foundation for investigating possible causes of glacier fluctuations in New Zealand. Here, I focus on moraines dating to the Lateglacial (15,000 to 11,500 years ago or 15 - 11.5 ka) and Holocene (last 11.5 ka). Second, the volume of present-day glacier mass balance and alpine climatological data in the Southern Alps is increasing. These data aid in improving numerical modelling simulations through model evaluation. Third, simple numerical models that take into account the most important mountain glacier processes now exist and improved computing power also makes it possible to carry out a large number of simulations. 2.1 Moraines and palaeoclimate In this thesis, I examine the often-held premise that glaciers are sensitive palaeoclimate indicators and that moraine sequences record past climate events. Moraines are landforms that have been deposited directly by glaciers, and they normally consist of unstratified sediment (Embleton and King, 1975). The form and composition of moraines reflect the local topography, lithology, climate, debris supply, and the efficiency of sediment transport from the glacier to the proglacial environment (Benn et al., 2003). Moraine size is related to glacier debris flux and the amount of time the glacier terminus remained at a single location (Evans, 2003). The term ‘moraine’ covers a wide range of depositional features, and moraines that formed at the snout of a glacier when it was stationary are ideal for chronologic dating and interpreting past climate (Embleton and King, 1975). 11 Interpretations of palaeoclimate from moraines have developed over the last 200 years. In 1837, Louis Agassiz was one of the first to propose that there had been a past ice age, and he observed moraine evidence from the European Alps, Scotland, and North America that supported his theory (Imbrie and Imbrie, 1986). Subsequent theories related summer insolation intensity to the pacing of the ice age cycles (Milankovitch, 1941; Hays et al., 1976), thus relating past climate to past glacier size. As a means of estimating past climate, glacial geologists have used constraints from glacial geomorphology, such as moraine location and other evidence where present, such as trimlines, to define past glacier extent. This allows equilibrium-line altitudes (ELAs) to be reconstructed using the accumulation-area ratio (AAR) method (e.g. Porter, 1975a; Broecker and Denton, 1990). In this method, the total area of the past glacier extent is divided into the accumulation and ablation areas (using a ratio of 2:3, accumulation:total area). The reconstructed ice surface elevation at the dividing line between the two areas is an estimate of the past ELA. The difference between past and present-day ELA (∆zELA ) can then be converted into an estimate of temperature change (∆T ) from present-day, , typically -6.5◦ C km−1 ): using a temperature lapse rate ( dT dz ∆T = ∆zELA dT dz (2.1) This approach is advantageous because of its simplicity, resource and time efficiency, and ability to provide reasonable ELA and temperature change (∆T ) estimates. Assumptions implicit in the AAR method, however, translate into considerable uncertainty in the ELA and ∆T results. In recent years, as moraine chronologies have improved in precision and coverage, and as present-day glacier-climate interactions are observed, some researchers have questioned the robustness of climate interpretations that can be made from moraines (Balco, 2009; Kirkbride and Winkler, 2012). For example, several modelling studies have shown that significant (>1 km) changes in glacier length can occur in response to interannual stochastic variability in temperature in an otherwise steady climate (Roe, 2011). This finding has 12 caused researchers to question whether or not moraines located close to the glacier terminus require a systematic change in climate in order to be deposited. To evaluate the robustness of palaeoclimate-moraine interpretations, these alternative ways of explaining changes in glacier length should be explored. The New Zealand landscape contains a variety of glacial deposits, including large, well-preserved moraines. Detailed geomorphic maps and moraine chronologies help to clarify the landscape history and allow for a comparison between moraine records and other palaeoclimate records. 2.2 New Zealand Lateglacial and Holocene palaeoclimate In general, data from terrestrial climate proxy records are more sparse in the Southern Hemisphere than in the Northern Hemisphere due, in part, to there being a smaller land area. New Zealand is one of the few locations in the southern middle latitudes with active mountain glaciers, which can provide a basis for understanding past glaciations. The palaeoclimate of New Zealand since the last glacial maximum (LGM) has been reconstructed from deposits left by glaciers, stable isotopes, flora, and fauna (e.g. McGlone, 1995; Fitzsimons, 1997; Singer et al., 1998; Newnham and Lowe, 2000; Turney et al., 2003; Williams et al., 2005; Anderson and Mackintosh, 2006; Hajdas et al., 2006; Alloway et al., 2007; Vandergoes et al., 2008; Williams et al., 2010). Some moraine sequences show a Lateglacial cooling associated with the Antarctic Cold Reversal (ACR), as well as several Holocene glacial advances. These events differ in timing and magnitude to those which occurred in Europe since the LGM (Gellatly et al., 1988; Ivy-Ochs et al., 1999; Schaefer et al., 2009; Putnam et al., 2012). The Lateglacial chron includes two abrupt (∼2 ka) cooling events. The Younger Dryas (YD) (12.7 - 11 ka) is unambiguously registered in the North Atlantic region in ice cores, pollen records, and moraines (EPICA Community Members, 2006; Broecker et al., 2010; Carlson, 2010). Another event, known as the ACR (15 - 13 ka), is primarily found in Antarctic ice cores, Southern Ocean marine sediment cores, and South American and New Zealand moraines and pollen records (Blunier et al., 1997; Newnham and Lowe, 13 2000; Fogwill and Kubik, 2005; Hajdas et al., 2006; Barrows et al., 2007a; Carter et al., 2008; Vandergoes et al., 2008; Kaplan et al., 2010; Putnam et al., 2010a). Learning more about the YD and ACR climate signals could help to clarify interhemispheric heat transfer, abrupt climate drivers, and responses of the ocean, atmosphere, and terrestrial systems to climate change. The timing of a Lateglacial ice advance in the central Southern Alps was originally thought to be synchronous with the YD found in North Atlantic / European records (Denton and Hendy, 1994; Ivy-Ochs et al., 1999) but has now been shown to be asynchronous with this event (Turney et al., 2007; Kaplan et al., 2010; Putnam et al., 2010a). Studies of glaciers and their moraine records, pollen, and chironomids in the Southern Alps suggest that a cooling occurred during the ACR (Newnham and Lowe, 2000, 2003; Vandergoes et al., 2005, 2008). After this cool period, glaciers in New Zealand retreated during the YD, while glaciers in Europe readvanced (Golledge et al., 2009; Kaplan et al., 2010). A previous attempt to model New Zealand glacier extent during the Lateglacial was carried out for Franz Josef Glacier (Anderson and Mackintosh, 2006). The Waiho Loop moraine, which has been associated with the Lateglacial (Denton and Hendy, 1994; Turney et al., 2003), represents a 10 km advance of the Franz Josef Glacier from its modern position (Figure 2.1). Anderson and Mackintosh (2006) discovered that a 4.1 - 4.7◦ C drop in mean-annual temperature, 400% increase in mean annual precipitation, or some combination of the two would be necessary to explain this advance (Anderson and Mackintosh, 2006). However, the origin (Tovar et al., 2008) and age (Barrows et al., 2007b) of this moraine have been questioned, causing some doubt about the palaeoclimate interpretations that can be drawn at this location (Evans, 2008). Unlike those in the European Alps, glaciers in New Zealand appear to have retreated steadily during the Holocene, making the Little Ice Age a relatively minor Holocene glacial event (Gellatly, 1984; Gellatly et al., 1988; Ivy-Ochs et al., 1999; Schaefer et al., 2009; Putnam et al., 2012). As a result, New Zealand glaciers did not override and destroy earlier Holocene moraines. This asynchronous behaviour between hemispheres has been attributed to climate processes associated with the bipolar seesaw, insolation, and 14 regional climate feedbacks (Broecker, 1998; Schaefer et al., 2009; Kaplan et al., 2010; Putnam et al., 2010a, 2012). Several modelling studies have shown that stochastic, interannual temperature variability can cause kilometre-scale glacier terminus fluctuations (Oerlemans, 2000; Roe and O’Neal, 2009; Roe, 2011), complicating the traditional view that moraines are deposited as a direct consequence of a change in climate. These simulations shed doubt on the need to explain Holocene glacier fluctuations by invoking a significant shift in climate (∼30year mean). In other words, these authors contend that glacier fluctuations of Holocene scale might reflect nothing more than interannual stochastic variability in climate. Consequently, the palaeoclimate significance of recently published New Zealand Holocene moraine chronologies (Schaefer et al., 2009) have been questioned (Balco, 2009). In summary, the New Zealand landscape contains many moraine sequences that can be mapped, dated, and potentially used to interpret past climate, along with other climate proxy records. There are, however, uncertainties and assumptions made in extracting palaeoclimate data from moraine records, especially when applying simple glacier lengthtemperature change relationships such as the AAR method. The relationships between glacier length and climate change should be further examined using models of glaciers to test these assumptions and to better understand what moraine sequences can tell us about past climate. 2.3 Lateglacial and Holocene glacial geology of New Zealand The extensive nature of many New Zealand moraines is in part due to the high sediment flux delivered to the glacier surface. The central Southern Alps of New Zealand provide an ideal environment for generating and supplying ample debris to glaciers (Shulmeister et al., 2010b). This is due to the combination of steep valley walls, substantial exhumation rates related to the vertical velocity of surface uplift, a highly-fractured and low-grade metamorphic lithology (‘greywacke’), large amounts of orographic precipitation, and associated high denudation rates (Little et al., 2005; Houlié and Stern, 2012). The high 15 debris supply, along with other conditions mentioned above (Section 2.1), typically result in voluminous moraines (Benn et al., 2003). Lateglacial and Holocene moraines are located between the LGM and present-day ice extents (LGM marked by a pink line and light blue shows present-day ice extent in Figure 2.1). The best-preserved moraines in New Zealand have been the focus of decades of mapping and dating research, and have been used to evaluate global-scale climate events (Suggate, 1965; Soons and Gullentops, 1973; Burrows, 1975; Mercer, 1984; Gellatly et al., 1988; Fitzsimons, 1997; Ivy-Ochs et al., 1999; Schaefer et al., 2009; Kaplan et al., 2010; Shulmeister et al., 2010a; Putnam et al., 2012). 170°0'0"E 43°0'0"S ¯ 171°0'0"E 43°0'0"S 3 2 4 1 7 6 5 44°0'0"S Rakaia River Lake Tekapo 44°0'0"S Lake Pukaki Lake Ohau 0 5 10 15 km 170°0'0"E LGM ice extent Rivers Lakes Glaciers 171°0'0"E Figure 2.1 Map of shaded topography of central South Island, New Zealand with the last glacial maximum ice extent (pink outline Barrell et al. (2011)) and the Main Divide of the Southern Alps (dashed red line). Numbers mark locations of 1) Franz Josef Glacier, 2) Hokitika, 3) Arthurs Pass, 4) Christchurch, 5) Irishman Stream, 6) Whale Stream, and 7) Cameron valley. 2.3.1 Mapping Lateglacial and Holocene moraines Until recently, glacial deposits in different regions of the South Island (Mackenzie basin, West Coast, Rakaia, etc.) were mapped and described by different geologists and then 16 correlated (Gage and Suggate, 1958; Suggate, 1961, 1965; Gair, 1967; Soons and Gullentops, 1973; Porter, 1975a; Shulmeister et al., 2001; Suggate and Almond, 2005; Barrell et al., 2011). Moraines from different areas were matched using the degree of preservation, geographic position, and relative weathering of boulders (e.g. quartz vein heights, weathering rinds, and lichenometry). More recently, the central part of the South Island’s glacial geomorphology has been systematically mapped and presented in a monograph by GNS Science (Barrell et al., 2011). Glacial deposits in these maps are categorised by relative age, with colour used to separate Holocene, Lateglacial, LGM, and pre-LGM deposits. The maps are digitised and provide larger-scale context to guide and focus the SED investigations, which in turn, provide absolute age control. Higher resolution segments of the central South Island glacial geomorphology map have accompanied recent publications, detailing the cosmogenic moraine chronologies by the Denton group (Schaefer et al., 2009; Kaplan et al., 2010; Putnam et al., 2010a, 2012, 2013). The maps provide the context of each moraine within a sequence, showing moraine ridge continuity, size, and complexity. These moraine maps are ideal for modelling studies, where glacier simulations are compared to mapped moraine positions. 2.3.2 Dating of moraines Moraine dating, when applied at a global scale, allows for comparisons of climate events, which may reveal large-scale climate teleconnections (Clark et al., 2009). Determining the age of a moraine is the first step toward the goal of assessing its regional or perhaps, global significance. Initial dating surveys used weathering rind thicknesses (Gellatly, 1984), lichenometry, and Schmidt hammer tests (Winkler, 2005) to compare deposits from valley to valley. However, these techniques provide calibrated ages, rather than absolute dates and their associated uncertainties in age do not allow for broader-scale comparisons. The first radiometric ages of moraines in New Zealand, and elsewhere, started with radiocarbon dating of buried organic material (Burrows, 1988; Lowell et al., 1995). This method works well in areas with high amounts of organic material (i.e. forested areas). 17 For example, radiocarbon ages from tree bark buried under till on a bedrock knoll on the West Coast were interpreted to represent the timing of an advance of Franz Josef Glacier to the Waiho Loop moraine (Denton and Hendy, 1994). More recently, cosmogenic SED has provided the means to accrue a large number of moraine ages in the drier, organicscarce regions east of the Main Divide of the Southern Alps (henceforth referred to as the main divide, red dashed line in Figure 2.1). Improvements in SED techniques (i.e. refined boulder sampling, low-background carrier, and high accelerator mass spectrometer currents) have allowed for more precise moraine chronologies (Schaefer et al., 2009), extending over a timescale from millions of years right up to the last century. In addition, a local 10 Be calibration site in the Macaulay valley (Putnam et al., 2010b) has improved dating accuracy in New Zealand. Sampling campaigns in the central South Island targeted areas with multiple moraine ridges and little post-depositional disturbance. Such areas are readily identified using aerial photographs and detailed glacial geomorphology maps (Barrell et al., 2011). Emerging chronologies of the Lateglacial and Holocene moraines in the Southern Alps help to fill gaps in knowledge of Southern Hemisphere glacier fluctuations (Schaefer et al., 2009; Kaplan et al., 2010; Putnam et al., 2010a; Kaplan et al., submitted; Putnam et al., 2012), dramatically improving our understanding of past glacier behaviour in New Zealand. Thus far, the published chronologies from the Denton group have highlighted moraine ages, which are interpreted as representing the time when the glacier pulled away from its steady-state position at the moraine, dating from the penultimate glaciation (Putnam et al., 2013) to the late Holocene (Schaefer et al., 2009; Putnam et al., 2012). Moraines dated to the Lateglacial suggest that glaciers advanced coincident with the ACR (∼13 kya) and subsequently retreated during the YD (Kaplan et al., 2010; Putnam et al., 2010a). Holocene moraine sequences show multiple glacial advances (e.g. ∼8,200, ∼6,000, and 500 years ago) with an overall retreat to the present-day ice extent (Schaefer et al., 2009; Putnam et al., 2012). 18 2.4 Climate of the Southern Alps New Zealand is positioned in the southwest Pacific, just north of the present-day location of the Sub-Tropical Front (Carter et al., 2004). The mountains of the South Island trend northeast along the Australia-Pacific plate boundary zone, with many peaks reaching between 2000 and 3000 m asl in the central Southern Alps (highest peak, Aoraki / Mt. Cook 3754 m asl). These mountains intersect the southern middle latitude westerly winds, creating a strong orographic effect and precipitation gradient from west to east (Figure 2.2). The pattern of mean annual precipitation is imprecisely known, but is thought to peak at ∼10 m (Griffiths and McSaveney, 1983; Henderson and Thompson, 1999) between the Alpine Fault and main divide and decrease almost exponentially with distance east over to the Canterbury plains (Chinn and Whitehouse, 1980; Henderson and Thompson, 1999; Salinger and Mullan, 1999). The maritime climate is characterised by a low seasonal range in temperatures, which means glaciers with high elevation catchments and low elevation tongues can receive accumulation as well as ablate ice year-round. 2.4.1 Southern Alps climate variability Synoptic-scale atmospheric and oceanic circulation in the South Pacific region ultimately determine the present-day climate variability in New Zealand (Kidson, 1994; Chinn, 1995; Mullan, 1995; Salinger and Mullan, 1999; Clare et al., 2002). The synoptic circulation is responding to larger-scale climate patterns and modes, such as the El NiñoSouthern Oscillation (ENSO), Southern Annular Mode (SAM), and Interdecadal Pacific Oscillation (IPO). These modes or oscillations all affect the precipitation and temperature variability observed in New Zealand to some degree (Jiang et al., 2012). Climate data from Hokitika (West Coast), Arthurs Pass (main divide), and Christchurch (east coast) provide an idea of present-day interannual precipitation and temperature variability (Table 2.1). Proxy records from New Zealand indicate that synoptic circulation was also important during the Holocene (Ackerley et al., 2011; Lorrey et al., 2012; Fowler et al., 2012), but 19 170°E 171°E Precipitation (mm / a) 0 - 2,000 2,000 - 4,000 4,000 - 6,000 6,000 - 8,000 8,000 - 11,000 44°S 44°S ¯ 0 5 10 170°E 20 km 171°E Figure 2.2 Map of 30-year mean annual precipitation surface (Stuart, 2011), overlain on shaded topography of central South Island, New Zealand and the Main Divide of the Southern Alps (dashed black line). In the central Southern Alps, most of the annual precipitation falls west of the main divide due to a strong orographic effect. 20 there is still much left to learn. The frequency or intensity of the different climate modes were most likely different from today. For example, the onset of ENSO is thought to have occurred during the mid-Holocene, likely due to a change in orbital configuration. However, the frequency of El Niño / La Niña events was not the same as today (Clement et al., 2000; Liu et al., 2000; Moy et al., 2002). Table 2.1 Climate information from Hokitika, Arthurs Pass, and Christchurch stations. Statistics from daily mean precipitation and temperature data from 1980 to 2011 (NIWA, Retrieved 2009-2011). Hokitika (west coast) and Christchurch (east coast) data are presented for near sea-level west to east comparisons. Elevation Precipitation Mean (mm) Annual std dev. (mm) Annual std dev. (%) Temperature Mean (◦ C) Annual std dev. (◦ C) Annual std dev. (%) Daily std dev. (◦ C) 2.4.2 Hokitika 39 m asl Arthurs Pass 738 m asl Christchurch 37 m asl 2886 329 11 4661 726 16 596 113 19 12 0.5 4.1 1.84 7.7 0.51 6.6 2.30 11.6 0.42 3.6 2.75 Alpine meteorology While decades of weather station data have helped to describe the spatial and seasonal patterns of temperature and precipitation near sea level in New Zealand, alpine meteorology has received less attention. The relatively few high-altitude automatic weather stations (AWSs) that do exist are monitored by organisations such as universities, Crown Research Institutes (e.g. GNS, NIWA), recreational providers (e.g. ski field), or hydroelectric companies (e.g. Meridian). The stations (at least 11 stations monitored by the National Institute of Water and Atmospheric Research (NIWA) alone) measure air temperature, relative humidity, wind speed and direction, air pressure, solar radiation, reflected radiation, and sometimes precipitation (Gillett and Cullen, 2011). Maintaining an alpine AWS on a glacier surface can be difficult and no permanent glacier AWSs exist: ‘Automatic recorders are expensive however, and need to be robust, kea proof (at least in 21 New Zealand) and capable of being retrieved from under perhaps 20 m of snow’ (Bishop and Forsyth, 1988). Temperature lapse rates, which describe the change in temperature with increasing elevation, are affected by processes in the atmospheric boundary layer (Barry, 1990, 1992) and vary in time and with location (Blandford et al., 2008). However, due to the scarcity of high altitude climate station data and uncertainty in lapse rate variability, a single value is often used, typically an ‘environmental’ adiabatic lapse rate (-6.5◦ C km−1 ), which represents conditions roughly mid-way between the wet and dry adiabatic lapse rates (Sturman and Tapper, 1996). Lapse rates in New Zealand can be estimated using an empirical relationship between the lapse rate and the site latitude, elevation, and distance to the nearest coast (Norton, 1985). However, using a lapse rate from site-specific AWS measurements is more desirable if such data exist. 2.5 Glacier monitoring Glaciological studies, together with studies of alpine meteorology, help to establish links between characteristic weather patterns and glacier mass balance (Oerlemans, 1992; Oerlemans and Knap, 1998; Klok and Oerlemans, 2002; Hock, 2005). Glaciers fluctuate in size to equilibrate with their climate where ice extent is dependent on a combination of factors, including; ice flow rates, topography, total annual accumulation of snow, and total annual melt (Paterson, 1994). Although annual net mass balance (accumulation minus melt) varies from year to year, ice extent adjusts over longer timescales (years to centuries), thus glacier length changes typically mimic a smoothed climate signal. To interpret past climate from dated moraine positions requires an understanding of how glaciers in New Zealand respond to modern climate. Specifically, observations of glaciers in the Southern Alps provide a basis for tuning and evaluating glacier models. Glacier mass-balance studies in New Zealand have focused on several glaciers chosen for their accessibility, size, and importance. Long-term investigations have focused primarily on the Ivory, Tasman, Franz Josef, Brewster, and Fox glaciers (Anderton and Chinn, 1978; 22 Anderson et al., 2008, 2010; Gillett and Cullen, 2011; Purdie et al., 2011b). Brewster Glacier is currently a major part of the mass balance programme, with data reported annually to the World Glacier Monitoring Service (WGMS, 2011). One- and two-year mass-balance studies have been carried out on other South Island glaciers, including the Park Pass, Glenmary (Stumm, 2011), and Dart Glacier (Bishop and Forsyth, 1988). The Rolleston Glacier mass-balance programme started in November 2010 and regular monitoring is planned to continue into the future (personal communication T. Kerr, 2010). The end of summer snowline (EOSS) survey is the longest record of glacier observations in New Zealand (Chinn, 1995; Chinn et al., 2006, 2012). This survey includes 50 index glaciers from Fiordland to the Kaikoura Ranges and has been on-going since 1977 (Chinn, 1995; Willsman et al., 2010). The record is comprised of oblique aerial photographs, taken from roughly the same elevation and vantage point each year, that show the distribution of snow at the end of the mass balance year. The EOSS in the photograph is then visually compared to a topographic map to estimate the annual ELA. The EOSS at any one time varies geographically due to accumulation and ablation of snow (Chinn, 1995) ranging in elevation from ∼1600 m asl in the west to 2200 m asl in the east. Annual EOSS elevations of individual glaciers are highly correlated to the annual mean EOSS elevation of the Southern Alps, demonstrating that they largely respond to the same climate signal (Chinn et al., 2006, 2012). Several of the large, low-angle glaciers east of the main divide now show significant debris cover (Reznichenko et al., 2010; Anderson and Mackintosh, 2012), and many are calving into proglacial lakes that have formed as a consequence of recent thinning (Dykes and Brook, 2010; Chinn et al., 2012). Thick debris cover can insulate the ice below, decreasing the melt rates, but does not stop melt from occurring (Anderson and Mackintosh, 2012; Gardelle et al., 2012), whereas iceberg calving into lakes can increase the rate of glacier retreat. Both of these processes are imperfectly understood in the present-day and the impact of these processes is largely unknown for the past. Modelling sites were chosen to avoid glaciers with these features and associated, more complicated, responses to climate. Overall, New Zealand glaciers can be characterised as temperate mountain glaciers in a 23 wet, maritime environment that are highly sensitive to changes in climate due to their high accumulation and high ablations rates (high mass balance gradients). High surface angle glaciers, often but not exclusively found on the West Coast, have faster reaction times to climate changes than the large, low-angle glaciers typically found east of the main divide (Oerlemans, 2000; Chinn et al., 2012). 2.6 Glacier modelling Models are used to represent a complex natural system in a simplified manner. The most important components of a system are included in equations or parameterisations that attempt to describe the physical relationships. Glacier models vary widely in complexity to suit the requirements of the problem being investigated. Glacial geologists typically use the AAR method as a simple model, which can be applied to any catchment with a well-defined terminal moraine using a topographic map. Estimates of past ELA change can be made by outlining past glacier area, as determined by the moraine, and using an assumed AAR (typically of 0.6 to 0.66), where the ELA is solved for as the dividing line between the accumulation and ablation zones. Porter (1975a) used the AAR method to work out the ELA lowering during the Lateglacial period in New Zealand. An ELA lowering of 500±50 m was calculated for the Birch Hill moraine in the Tasman-Pukaki basin. The difference in past and present-day ELAs can then be converted into temperature change by using a temperature lapse rate (typically -6 to -6.5◦ C km−1 ). This method is fast and gives an estimate of temperature change for a steady state ‘snapshot’, but cannot account for mass balance influences such as topographic shading, hypsometry, accumulation concentration, nor distortions resulting from differences in glacier flow. Due to the many assumptions made, the results from this type of model include errors that are difficult to quantify (Plummer and Phillips, 2003; Oerlemans, 2005). Numerical models may offer more comprehensive palaeoclimate reconstructions than an AAR approach (Hubbard, 1997; Plummer and Phillips, 2003). Past glacier extents have 24 been successfully reproduced in a variety of locations, using models of intermediate complexity (Oerlemans, 1988; Mackintosh et al., 2002; Kessler et al., 2006; Golledge et al., 2009). These models typically consist of a mass balance component and an ice-flow component, each containing a suite of tunable parameters. In Chapter 3, I describe several different model types and the benefit and assumptions associated with each. To date, published palaeoclimate estimates from numerical simulations of past glaciers in New Zealand have come from a range of model types. These types include a degree-day model with a 1-D flowline of Franz Josef Glacier (Anderson and Mackintosh, 2006), a 3-D shallow-ice approximation (SIA)/shallow shelf approximation model for the entire Southern Alps (Golledge et al., 2012), a 2-D SIA simulation of the Rakaia catchment (Rowan et al., in press), and flowline models used to study LGM ice extent in the Pukaki and Ohau catchments (McKinnon et al., 2012; Putnam et al., 2013). Herman and Braun (2008) also produced a simulation of glaciation of the Southern Alps, however, their aim was to simulate landscape evolution. 2.7 Research questions Moraine records have been used as proxy indicators of past climate for more than a century. Models provide a potentially useful means to add value to the mapping and dating of glacial landforms. The following research questions and objectives stem from the desire to understand what moraines and their chronologies represent in terms of palaeoclimate and glacier response. The first two questions are related to past climate, the next three explore some finer details of glacier mass balance variability, and the last three concern the methodology employed. 1. What magnitude of past temperature and precipitation change do Lateglacial- and Holocene-aged moraines in New Zealand represent? 2. Are these estimates of past temperature and precipitation coherent with other climate proxy records in New Zealand? 25 3. To what extent does the re-mobilisation of snow by avalanches affect the presentday mass balance of New Zealand glaciers? 4. Does including a model of snow redistribution in a glacier model change the palaeoclimate reconstructions made? 5. Could natural, stochastic variability in climate have caused the Holocene glacier fluctuations? 6. What parameters are the modelled glaciers most sensitive to, and how does parameter choice affect the climate reconstructions? 7. How do the palaeoclimate estimates derived from modelling glaciers differ from previous estimates using simpler methods? 8. What can we learn from dynamic glacier models that we could not gain from using simpler methods? The overall aim of my work is to improve palaeoclimate estimates from Holocene and Lateglacial moraines in New Zealand. Additionally, I hope to provide an idea of how quickly and by what magnitude these glaciers can respond to changes in climate. Understanding the sensitivity of these glaciers allows for more insightful interpretations of moraine sequences and ultimately, climate. 26 Chapter 3 Methodology “Essentially, all models are wrong, but some are useful” Box and Draper (1987). The study of glaciers and their response to climate change has progressed greatly in the past few decades. Models have been an essential part of developing the understanding of glacier-climate relationships for the present-day, where glaciological and meteorological data are available. Numerical models provide the opportunity to better understand glacier-climate relationships and test different theories about what caused past glacier fluctuations. Models range in complexity, but they all attempt to represent the most important processes within a system, with a series of equations. Simplifying any natural system into a model results in an incomplete representation of the system, and the associated model assumptions and uncertainties must be discussed and explored (Oreskes et al., 1994; van der Veen, 1999). Numerical models that describe the physics of ice-flow and glacier mass balance can be used to simulate past glacier extents, as delimited by moraines. Such simulations improve our understanding of past climate, glacier-climate feedbacks, glacier response times, and ice thickness estimates, even in areas that lack geomorphic evidence. Simpler models, such as the accumulation-area ratio (AAR), cannot offer these insights, which is why I used numerical models to infer climate-glacier-moraine relationships. Here, I describe a range of available glacier models and explain why I chose the ones used in this study. I also discuss the limitations and advantages of a numerical modelling approach and explore the parameter values used in each model and in different chapters (Tables 3.2, 3.1, and 3.3). 27 3.1 Modelling approaches There are several reasons why I chose a numerical modelling approach to evaluate past climate from glacier records. First, input climate data and parameters are adjustable, which allows for testing of different palaeoclimate scenarios, such as precipitation change from present-day, snow avalanche accumulation, or natural stochastic variability in climate. Second, the adjustable parameters allow for sensitivity tests, such as assessing the glacier length fluctuations that result from step changes in temperature, that enhance the understanding of glacier-climate relationships. Third, the models that I used consider locationspecific features, such as aspect, bedrock slope, surface elevation, and topographic shading. Finally, the model can produce transient simulations, thereby allowing the evolving glacier response to a changing climate to be assessed. The three model components used in this thesis describe (1) surface energy balance, (2) gravitational snow mass transport and deposition (MTD), and (3) ice-flow. The distributed energy-balance model (EBM) accounts for shading, changes in albedo, and aspect, among other factors that affect melt. The gravitational snow MTD model uses the grid cell slope values and total snowfall to simulate small, frequent snow avalanches. The grid of redistributed snow is used to improve mass balance estimates. The ice-flow model calculates the 2-dimensional flow of ice across a landscape, which allows the size of the glacier to evolve as the mass balance changes. These models can be efficiently applied to small areas (<20 km2 ) with high-resolution digital elevation models (DEMs) (25 to 100 m). The mass balance distribution is ultimately what drives the glacier to advance or retreat. The close coupling between the models allows the responses and feedbacks between mass balance and evolving ice geometry to be explicitly considered. This coupling allows for the EBM to update accumulation and ablation rates for a particular surface topography (bedrock + ice thickness) because mass balance values in a cell are affected by changes in aspect, surface elevation, slope, or total ice cover. The ice thickness is also updated with the calculated horizontal ice flux from cell to cell at each timestep. Details about model setup, why I chose to use certain parameter values, and how these 28 models were applied to answer the research questions (Section 2.7), are explained in the following sections. All models use SI units, however some of the values in the following tables show non-standard units for easier comprehension. 3.2 Energy-balance model A variety of climate-glacier models exist, each with its own level of complexity and parameterisations. Positive degree-day (PDD) models are based on an empirical relationship between ice ablation and air-temperatures (Ohmura, 2001). A PDD model simplifies complex processes that are described in an EBM, and calculates melt by multiplying the sum of PDD with a PDD factor. This attempts to account for the surface energy balance processes that affect melt, including the incoming shortwave and longwave radiation and the turbulent heat fluxes (Braithwaite, 1995; Blard et al., 2011). The EBM explicitly considers the influences of topography on the surface energy balance within the atmospheric boundary layer (Equation 3.1). A surface EBM was chosen over a PDD model for several reasons. First, the EBM allows me to test the sensitivity of glacier extent to climatic inputs other than temperature. One aim of this thesis is to examine the relationship between glaciers and temperature change and test whether past glacier extents in New Zealand can be explained by other climate variables. To do this, I require a model that is not solely dependent on temperature. Second, although the EBM requires more input data (i.e. wind velocity, relative humidity, cloud cover, debris cover), which are poorly known in the past, I assume these data to be the same between past and present. Differences in these values (absolute or spatial gradients) would be interesting to investigate, but is beyond the scope of this study. Third, the EBM used here could be adapted to include the MTD model, which improves the snow distribution and mass balance gradient of an avalanche-fed glacier. The EBM uses a surface elevation grid (digital elevation model plus an ice thickness grid) and daily climate data. The climate data came from the National Institute of Water and Atmospheric Research (NIWA) CliFlo Database, NCEP reanalysis, and an annual mean precipitation surface (Stuart, 2011), empirically derived from the most comprehensive 29 network of precipitation stations available. Precipitation amount depends on the daily or monthly mean interpolated precipitation grids within my model and the annual mean precipitation surface (Stuart, 2011). Snowfall occurs in a grid cell in the model when the atmospheric temperature (Ta ) at that grid cell is at or below a snow/rain temperature threshold (Tsnow ). Tsnow was set to 1◦ C (274.15 K) for all simulations (Anderson et al., 2006) (Table 3.1) and varied in the sensitivity tests by ±1◦ C. If simulations include the MTD model, modelled snow is transported from steep slopes to lower, gentler slopes (see Section 3.3). Temperature lapse rates in New Zealand vary spatially and temporally, however, there are limited temperature measurements from alpine sites. An automatic weather station was installed in Irishman basin for a separate study (K. Sattler and others, Victoria University of Wellington) and I calculated monthly mean temperature lapse rates from March 2010 to January 2011 temperature data. There is a noticeable seasonal cycle in the lapse rate values (Figure 3.1a), with lower values in winter (July -3.5◦ C km−1 ) and higher values in summer (January -6◦ C km−1 ). After publication of this article, we had access to the 20112012 temperature measurements which showed a similar pattern in temperature lapse rate variability (Figure 3.1b). The monthly lapse rates calculated for 2010 to 2011 seem to work well for the 2011 to 2012 year (Figure 3.1b, red line) however, we offer an improved fit for the final three months using values of -6.5, -6.5, and -7.5◦ C km−1 for October, November and December, respectively. A DEM (GeographiX (NZ) Ltd) derived from elevational contours by Land Information New Zealand (NZMS260), provided surface elevation for all model simulations. The gridded horizontal resolution of the DEM applied differs between chapters, depending on the research question and domain size (Table 3.4). In Chapter 6, the DEM was altered by removing the present-day ice thickness, which was estimated using the ice surface slope and modelled mass balance. The estimation of Cameron Glacier ice thickness was also guided by ground-penetrating radar (GPR) data which showed the bed topography along the glacier’s central flowline (Appendix B). The gridded topography was then smoothed using an average block size of five grid cells (250 m) (Chapter 6). The ice-removal and topographic smoothing is important for minimising the ice flux correction in the ice-flow 30 8 Temperature (deg C) 6 8 a) 6 4 4 2 2 0 0 -2 -2 -4 -4 -6 -8 Mar Measured Temp Modelled Temp (variable lapse rate) Modelled Temp (lapse rate -6°C km-1) Apr May Jun Jul Aug Sep Oct Nov Dec Jan b) Measured Temp Modelled Temp (var lps rate) Modelled Temp (lps -6°C km-1) Modelled Temp (new var lps) -6 -8 Feb Mar Month of the Year (2010-2011) Apr May Jun Jul Aug Sep Oct Nov Dec Jan Month of the Year (2011-2012) Figure 3.1 Modelled versus measured temperatures at the location of a climate station in Irishman basin. Modelled temperature values (daily interpolated temperature, T , smoothed with a thirty-day running mean) were interpolated from climate stations at <800 m asl over the 2010-2011 (a) and 2011-2012 (b) intervals. Measured temperature data (black line) are compared to modelled temperatures that were calculated using the 2010-2011 monthly varying lapse rate (red line), using a new variable lapse rate (purple line, panel b only) and a constant lapse rate (blue dashed). Using a constant lapse rate achieves a poor match, especially during winter months (Jun-Sep). model (Section 3.4.1). Debris cover was mapped from aerial imagery taken on 30 January, 2006, provided by Terralink International. Modelled melt in cells with prescribed debris cover was reduced by ∼90%, based on previous estimates on New Zealand glaciers (Anderson and Mackintosh, 2012). A debris cover grid was only included in the present-day model run because the debris coverage is unknown for the past. It is likely that present-day debris cover is greater than in the recent past, due to glacier thinning (Stokes et al., 2007). The timestep of the EBM is adjusted depending on the question being addressed. For palaeo-glacier simulations, model runs were carried out for over >200 model years to ensure that the glacier reached equilibrium. The climatic input data for these simulations are monthly means based on present-day climate (1981-2010). The EBM uses a monthly timestep for these model runs, which allows the seasonal cycle to be captured. The EBM in Chapter 5 (snow transport simulations) uses daily climatic input data from various meteorologic stations over the 2009-2010 mass balance year that were interpolated across the 31 model domain. One goal of the daily simulations in Chapter 5 was to replicate measured snow depths with the model. This investigation required a daily timestep in the EBM to resolve snow deposition and melt during the 2009-2010 mass balance year. To simulate mass balance gradients seen today on the Cameron Glacier using a monthly timestep, the temperature input data were further adjusted in Chapter 6. The 30-year monthly mean (1981-2010) temperature input data were converted into ‘daily’ temperatures by adding them to a randomly generated number (normal distribution with a standard deviation of 2.56◦ C). The standard deviation of daily temperature variations in a month was extracted from the CliFlo daily temperature statistics (Mt. Potts Station, May 2009September 2012). Inclusion of this ‘daily’ temperature variability was important for capturing the turbulent heat fluxes with monthly-mean temperatures. The model calculates the energy available for ablation (QM ) using the following equation: QM = I(1 − α) + L ↓ +L ↑ +QH + QE + QR + QG (3.1) where I is the incoming shortwave radiation, α is the surface albedo, L ↓ is the incoming longwave radiation, L ↑ is the outgoing longwave radiation, QH and QE are the sensible and latent heat fluxes, QR is the heat supplied by rain, and QG is the ground heat flux. All heat exchanges are in units of W m−2 . Positive values of radiative and turbulent fluxes indicate energy available for melting snow and ice, and negative values indicate a loss of energy. QR was calculated assuming that precipitation temperature is equal to air temperature (Oerlemans, 1992) and QG was set at 1 W m−2 (Neale and Fitzharris, 1997). I will now discuss each of the main components of the energy balance equation (Equation 3.1). 3.2.1 Shortwave radiation Most of the heat available at the Earth’s surface is dependent on the Sun’s energy. Solar radiation lies almost entirely between 0.15 and 4 µm, whereas terrestrial radiation is 32 infrared and lies between 4 and 120 µm, hence separating them into shortwave and longwave radiation. The amount of shortwave radiation that influences a location is primarily dependent on the time of day, time of year, changes in Milankovitch cycles, latitude, topographic shading, and cloudiness. Incoming shortwave radiation (I) contains direct and diffuse components (Oerlemans, 1992) from calculated radiation at the top of the atmosphere (Eisenman and Huybers, 2006). (3.2) I = ta tc (Idif + Idir ) where ta and tc are transmissivity for air and clouds respectively, and Idif and Idir are the diffuse and direct components of insolation respectively (Oerlemans, 1992; Anderson et al., 2010). The transmissivity for air was calculated using:  ta = (0.79 + 0.000024z) 1 − 0.08 π2 − ϕsun π 2  (3.3) where z is the surface elevation (m asl), and ϕsun is the solar azimuth angle (radians). The transmissivity for clouds was calculated using: tc = 1 − (0.41 − 0.000065z)n − 0.37n2 (3.4) where n is cloudiness. The cloudiness parameterisation follows that used in Anderson et al. (2010) where a fitted cubic relationship describes the relationship between cloudiness and the fraction of clear-sky to measured radiation (Hock, 2005). For clear-sky conditions Idif is typically ∼15% of the total incoming shortwave radiation and for overcast conditions ∼85%. Diffuse radiation was calculated using: Idif = [0.8 − 0.65(1 − n)] S sin 33 π 2 −Z  (3.5) where S is the solar constant and Z is the solar zenith angle. S is described as:  S = Ψa 1 + 0.034 cos  2πN 365  (3.6) where Ψa is the atmospheric clear-sky transmissivity and N is the day number (Oerlemans, 1992). For unshaded areas, direct radiation and the cosine of the angle of incidence between the slope normal and the solar beam (θ) can be calculated using: Idir = (0.2 + 0.65) ∗ (1 − n)S cosθ (3.7) cosθ = cos(β)cosZ + sin(β)sinZ cos(ϕsun − ϕslope ) (3.8) where β is slope and ϕslope is the slope azimuth angle (aspect). The amount of outgoing shortwave radiation depends on the surface albedo. The parameterisations used to estimate albedo are described below. 3.2.2 Albedo Two different albedo parameterisations were used in this thesis because two different timesteps are used. Simulations with a daily timestep (i.e. Chapter 5) used a mode where albedo decays with time (days) since the last snowfall (Oerlemans and Knap, 1998) (Table 3.4). This time-dependent parameterisation has been evaluated using measurements of incoming and reflected radiation data from Brewster Glacier, New Zealand (Anderson et al., 2010). All other simulations used a mode where albedo is dependent on snow thickness and the elevational difference from the modelled equilibrium-line altitude (ELA) (Oerlemans, 1992). This second mode is applied to models with monthly timesteps, where a time-dependent albedo would not be appropriate. The time-dependent albedo parameterisation depends on the snow cover and different 34 albedo values for fresh snow (αf rsnow ), firn (αf irn ), and ice (αice ) (Oerlemans and Knap, 1998) (Table 3.1). s (3.9) αsnow = αf irn + (αf rsnow − αf irn )e− t∗ where s is the time since the last snowfall occurred (days) and t∗ is a timescale (days) determining how fast the snow albedo approaches the firn albedo after snowfall (Oerlemans and Knap, 1998). αsnow is then used to calculate α: α = αsnow + (αice − αsnow )e− Dsnow d∗ (3.10) where Dsnow is snow depth (metres in snow-water equivalent or m s.w.e.) and d∗ is a characteristic scale for snow depth (m s.w.e.). A difference in albedo of redistributed snow in the avalanche model was not included. These schemes are all empirically derived, and the key parameters (e.g. αsnow ) are somewhat uncertain and therefore were evaluated in sensitivity tests. Alternatively, when using the ELA-dependent albedo parameterisation (Oerlemans, 1992), the ‘background’ albedo, αb , was first determined: 0.18 αb = 0.43 + arctan π  z − zELA + 300 200  (3.11) where z is surface elevation and zELA is the equilibrium line altitude (m asl). The albedo, α, is then calculated, where the albedo of snow, αsnow , is equal to 0.72 (Table 3.1). α = αsnow − (αsnow − αb )e−5Dsnow (3.12) For Chapters 4 and 6 and Appendix A areas of the glacier with no snow cover have an albedo set to the ice albedo value (αice , set to 0.34). 35 3.2.3 Longwave radiation Incoming longwave radiation comes from the surrounding terrain (e.g. heated rocky slopes) as well as the atmosphere (e.g. cloudiness and relative humidity). Unlike shortwave radiation, longwave radiation also depends on the air temperature (Equation 3.13). The viewfield for each cell specifies the amount of surrounding terrain and atmosphere at a particular location. The relative amounts are adjusted in response to ice thickness growth or decay. Incoming longwave radiation L ↓ is calculated using the Stefan-Boltzmann constant σ (J K−4 m−1 s−1 ), effective atmospheric emissivity ǫa , atmospheric temperature (approximated by Ta ), and the terrain viewfield Vf in the atmospheric (first) component, and, in addition, terrain emissivity ǫt and terrain temperature Tt in the terrain (second) component (Corripio, 2003; Plummer and Phillips, 2003; Anderson et al., 2010). L ↓= (σǫa Ta 4 )(Vf ) + (σǫt Tt 4 )(1 − Vf ) (3.13) Terrain temperature Tt is taken as atmospheric temperature for cells that are not snow- or ice-covered and as 273.15 K for cells with snow or ice. Terrain emissivity ǫt is set to 0.4 (Plummer and Phillips, 2003). Outgoing longwave radiation L ↑ is calculated using: L ↑= −ǫsnow σTsurf 4 (3.14) where ǫsnow is snow emissivity (set to 0.99) and Tsurf is the surface temperature (set to 273.15 K). I assume that the glacier surface is at the melting point, which is appropriate for temperature glaciers (Paterson, 1994; Ohmura, 2001), and this assumption makes the resulting L ↑ a constant equal to 317 W m−2 (Plummer and Phillips, 2003). 36 3.2.4 Turbulent heat fluxes Turbulent heat fluxes (QH and QE ), which can make up half or more of the energy available for melt in maritime environments (Anderson et al., 2010), were calculated using the bulk method (Oerlemans, 1992; Klok and Oerlemans, 2002; Oerlemans and Grisogono, 2002). QH = ρa cp kH Uw (Ta − Tsurf ) (3.15) where ρa is the density of ambient air, cp is the specific heat capacity of ambient air (J kg−1 K−1 ), kH is the bulk transfer coefficient for heat and Uw is wind speed (m s−1 ). QE = 0.622ρa kE Uw λv (q − qs ) p (3.16) where kE is the bulk transfer coefficient for vapour, λv is the latent heat of vaporisation for water (J kg−1 ), q is the vapour density of ambient air (kg m−3 ), qs is the vapour pressure on a melting surface (610.8 Pa), and p is the air pressure (Pa). k0 2    (1 − 5.2Rb )2  kH = z z log z0 log z0H (3.17) where k0 is the von Kármán’s constant (set to 0.4), z0 is the roughness length for wind (defined by Zsnow and Zice , m), z0H is the roughness length for sensible heat (m), and Rb is the Richardson stability criterion (Oerlemans, 1992; Anderson et al., 2010). The bulk transfer coefficient for heat (and similarly for vapour) depends on the roughness length for wind. Roughness lengths are difficult to measure and parameterise and have typically been used as a tuning parameter. Here we use previously published values for the roughness length of snow Zsnow and ice Zice (Table 3.1). The bulk transfer coefficients use the Richardson number stability criterion to keep the bulk Richardson number from being unrealistically high (Oke, 1987; Anderson et al., 2010). 37 In summary, the EBM used in this thesis was tailored to the specific research question in each chapter. The timestep was prescribed as daily or monthly for present-day or palaeosimulations, respectively. The albedo parameterisation and input climate data that were used in the model depended on the timestep. The ability of the EBM to accommodate different model resolutions, timesteps, and subroutines (such as the MTD model) allowed me to apply the EBM to all of the studies presented in this thesis. Table 3.1 Energy balance model parameters, values, and references. Parameter Snow/rain temperature threshold Temperature lapse rate Albedo ELA Roughness parameter for ice Roughness paramter for snow Albedo of snow Zice Zsnow αsnow Albedo of fresh snow Albedo of firn Albedo of ice Albedo characteristic depth Albedo characteristic time scale Maximum snow depth Model time step αf rsnow αf irn αice d∗ t∗ Dmax dt 3.3 Symbol Tsnow dT dz αELA Values 1◦ C variable ◦ C km−1 2100 m asl 0.004 m 0.001 m 0.72 0.75 0.9 0.53 0.34 11 mm s.w.e. 21.9 days 0.2 or 0.4 m s.w.e. variable Source (Anderson et al., 2006) refer to individual chapters (Oerlemans, 1992) (Anderson and Mackintosh, 2012) (Brock et al., 2006) (Oerlemans, 1992) Chapter 6 (Oerlemans and Knap, 1998) (Oerlemans and Knap, 1998) (Oerlemans and Knap, 1998) (Oerlemans and Knap, 1998) (Oerlemans and Knap, 1998) refer to individual chapters see Table 3.4 Gravitational snow mass transport and deposition model Snow avalanches contribute snow mass directly to glaciers located below areas of steep terrain (Lossev, 1967; Kotlyakov, 1973; Paterson, 1994; Purdie et al., 2011a). Small, frequent avalanches result in relatively thick snow deposits, typically near the base of steep valley walls, and can allow glaciers to persist even when their glacier surface is below the local ELA (Hewitt, 2005, 2011). Without avalanche processes, these low-elevation glaciers would not be sustainable and would retreat. Glaciers that receive avalanche accumulation are typically avoided for studying because of mass balance complexity as well as field hazards and difficult accessibility. However, avalanche-fed glaciers are common in the central Southern Alps where uplift, faulting, landsliding, and glacial and fluvial processes have exposed steep slopes that are prone to small, frequent snow avalanches. Because of this, the influence of avalanche accumulation on glaciers in New Zealand 38 needs to be better understood, especially in terms of how these glaciers response to climate change. Simple avalanche models have been used to simulate the effect of avalanche deposits on modelled glacier mass balance in Europe (Machguth et al., 2006; Dadic et al., 2008), and with proper tuning and evaluation, can be applied to studies in the Southern Alps. To parameterise small snow avalanches, I employ the gravitational MTD model following Gruber (2007). This parameterisation requires minimal input data, little computational time, and has been shown to produce appropriate snow distribution patterns on glaciers such as Vadret da Misaun (Gruber, 2007) and Claridenfirn (Machguth et al., 2006). The gravitational MTD used here was chosen over another snow avalanching model known as SnowSlide (Bernhardt and Schulz, 2010). This is because the Gruber (2007) model provides a larger number of adjustable and testable parameters, allowing for an improved fit between modelled and measured snow deposition. input slope limit: gamma: input: deposition limit: 75 50° 1 100 m swe 20 m swe 10 m swe 5 m swe 50 100 Elevation (m asl) Elevation (m asl) 100 25 75 50 a c 100 Distance (m) 150 slope limit: deposition limit: gamma: input: 75 50 200 50° 10 m swe 1 50 m swe 100 m swe 200 m swe 0 50 100 Elevation (m asl) 50 100 Elevation (m asl) 10 m swe 1 100 m swe 50° 40° 30° 25 0 25 100 Distance (m) 150 slope limit: deposition limit: input: gamma: 75 200 50° 10 m 100 m swe 1 0.5 0.1 50 25 b 0 deposition limit: gamma: input: slope limit: d 50 100 Distance (m) 150 200 0 50 100 Distance (m) 150 200 Figure 3.2 One-dimensional example of snow redeposition from the MTD parameterisation modified from Gruber (2007). Panels show avalanche profiles resulting from changing parameter values, a) increasing deposition limit (Dlim ) leads to thicker, less-extensive deposits, b) increasing total incoming snow will increase the extent of the avalanche deposit, c) increasing slope limit (βlim ) causes snow to stay close to the valley wall, and d) increasing gamma (γ) also results in less-extensive, thicker deposits. For these tests, the total incoming snow is a 100 m s.w.e. column deposited at the circle labelled ‘input’, which is an unrealistic process but helps to highlight the shape of the deposit. The MTD model is used in Chapters 5 and 6 and Appendix A. Chapter 5 explored param39 eter value choice to improve the fit between modelled and measured snow depths from a GPR survey. Chapter 6 and Appendix A use parameter values that were appropriate for the topography and model timestep. It is likely that many other glaciers in the Southern Alps of New Zealand are receiving accumulation through gravitational snow redistribution processes and this model could be used to improve the understanding of glacier mass balance and catchment-scale snow distribution patterns in other areas. I will now describe the MTD model components. The model begins with a subroutine designed to ‘fill’ any holes in the DEM to ensure that the topography is ‘drainable’. This process is also referred to as ‘iterative sink-filling’. There are two main, slope-dependent components in the gravitational MTD model; removal and deposition. The percentage of removed snow (or mobile mass Mm (m s.w.e.)) in the model is zero for cells with a slope of less than 40◦ and increases linearly to 100% for cells steeper than 70◦ . The Mm is the excess snow mass available for redeposition, and all remaining snow stays in place. The deposition mechanism computes the snow depth (Dsnow ) including the original snow cover and Mm using the specified Dmax curve: Dsnow    Mm =   D if Mm < Dmax (3.18) if Mm ≥ Dmax max where Dmax is the maximum amount of snow deposition allowed in a cell (m s.w.e.), which is ultimately slope-dependent (Equation 3.19). The simple function below describes the relationship between Dmax and the local slope angle β (◦ ): Table 3.2 Mass transport and deposition model parameter setup. Parameter Maximum snow depth limit (m s.w.e.) Maximum slope for snow deposition (◦ ) Slope influence 40 Symbol Dlim βlim γ Values 0.4, 0.6, 0.8 35, 40, 45 0.075, 0.1, 0.125 Dmax =   γ   (1 − β )Dlim βlim if β < βlim (3.19) if β ≥ βlim   0 where Dlim is the maximum deposition limit (where β= 0◦ ), γ alters the shape of the deposition-slope curve and βlim is the maximum slope at which snow begins to deposit. The parameters explored in this model alter the deposition-slope curve (Table 3.2) shown in Figure 3.3, where Dlim sets the y-axis intercept and βlim sets the x-axis intercept. The change in parameter values alters the shape and extent of snow avalanche deposits (Figure 3.2), which influences the goodness of fit between measured and modelled snow maximum deposition (m s.w.e.) depth and distribution (see Chapter 5). 0.8 0.7 0.8 a) Dlim Dlim = 0.4 0.7 Dlim = 0.6 Dlim = 0.8 0.6 0.8 b) g g = 0.075 g = 0.1 g = 0.125 Dlim = 0.4 b lim = 45 0.6 g = 0.125 b lim = 45 0.7 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 10 20 30 40 slope angle (°) 50 0 0 10 20 30 40 slope angle (°) b lim = 35 b lim = 40 b lim = 45 Dlim = 0.4 g = 0.125 0.6 0.5 0 c) b lim 50 0 0 10 20 30 40 50 slope angle (°) Figure 3.3 Plots of slope angle versus maximum deposition to show parameter influences on the slope-deposition curve. Parameters evaluated are a) Dlim (yellow circles) set to 0.4 (dotted line), 0.6 (dashed line), and 0.8 (solid line) m s.w.e., b) γ = 0.075 (dotted line), 0.1 (dashed line), and 0.125 (solid line), and c) βlim (green circles) set to 35 (dotted line), 40 (dashed line), and 45◦ (solid line). Without the MTD model, the snow accumulation in the EBM increases linearly with elevation as a function of the temperature lapse rate. Daily modelled snowfall depends also on precipitation amount, snow temperature threshold, and air temperature at that timestep. My GPR data show that simulated snow distribution on a glacier névé adjacent to steep slopes (>45◦ ) is greatly underestimated without the MTD model. The modelled end of summer mass balance for the Cameron Glacier, in the absence of the MTD model, shows ablation over the majority of the glacier surface (Figure 3.4c-d), which is consistent with inferences that this glacier lies largely below the local estimated ELA (2100 m asl) (Chinn, 1995; Willsman et al., 2010). In other words, without the avalanche-derived snow input, this glacier would not be able to sustain its current size. Modelling the snow depth 41 distribution appropriately is therefore critical for simulating the present-day glacier, as well as its past behaviour. b) a) Legend slope (° ) 0 - 30 30 - 40 40 - 50 50 - 60 60 - 70 50 m c) 0 0.150.3 d) 0.6 0.9 1.2 ¯ 1.5 km 10 5 0 -5 -10 -15 Annual mass balance (m for 2009-2010) 25 m -20 Figure 3.4 Cameron Glacier model domain showing a) slope (◦ ) at 25 m grid resolution, b) slope at 50 m grid resolution, c) modelled mass balance (m a−1 ) using the MTD model, and d) modelled mass balance without the MTD model for the 2009-2010 simulation. The black line represents part of the end of summer snowline measured on 19 March, 2010 using a hand-held GPS. I adjusted the Dlim and βlim parameter values from those used in Chapter 5 to values that enhance snow dispersion in Chapter 6. This adjustment helped to achieve a similar mass balance distribution between the EBM run at a monthly timestep, with a 50 m grid resolution, and the EBM run at a daily timestep, with a 25 m resolution, for the Cameron Glacier. When ice-flow dynamics are included in the monthly EBM, the ice-surface slope decreased through time (Chapter 6). Parameter adjustments also account for slightly lower slope values in the 50 m than in the 25 m grid (Figure 3.4a-b). Together these effects meant that the glacier was ‘flatter’ in a more subdued topography, and consequently, I changed the parameter values from βlim = 35◦ and Dlim = 0.4 m (Chapter 5) to βlim = 45◦ and Dlim = 0.2 m (Chapter 6) so that extensive avalanche deposits were still simulated. 42 3.4 Models that describe glacier flow Ice-flow models describe the mechanical movement of ice, based on the Navier-Stokes equation: ρice ∂U + ρice (U · ∇)U = −∇p + η∇2 U + F ∂t (3.20) is the rate of change of momentum, ρice (U · ∇)U is the convective force, where ρice ∂U ∂t −∇p is the pressure force, η∇2 U is the viscous force, and F is the external force. By neglecting the inertial terms (because of its slow movement) and assuming incompressability and a constant dynamic viscosity (µ), the Navier-Stokes equations for fluid flow reduces to the Stokes equations for creeping flow. The left-hand side of Equation 3.20 describes the inertial component and is neglected in the Stokes equations. 0 = −∇p + η∇2 U + F (3.21) The viscous force can be written in terms of the deviatoric stress tensor (τij′ ). Assuming incompressability (∇ · U = 0) and that the external force, F , primarily represents gravity (ρice g), the Stokes equation can be written as: 0 = −∇p + ∇ · τij′ + ρice g (3.22) To reduce computational time associated with the full-Stokes solution, many ice-flow models use a series of simplified equations, known as the shallow-ice approximation (SIA). The SIA was originally developed for describing large-scale flow dynamics in ice sheets, where the ice thickness is small compared to the horizontal extent (e.g. aspect ratio of glacier thickness to width is <0.2) (Leysinger Vieli and Gudmundsson, 2004) and the ice surface slopes and bed slopes are gentle (bed slope <11◦ ) (Le Meur et al., 2004; Egholm et al., 2011). The small aspect ratio is the dimensionless parameter upon which 43 the scale analysis is based. The applicability of this scaling depends on the similarity between the ratio of ice thickness to aerial extent and the ratio of vertical to horizontal velocities (Baral et al., 2001; Le Meur et al., 2004). This approximation is therefore less appropriate for ice divides or areas with steep bed slopes. In the most commonly used zeroth-order SIA, the influences of lateral drag and longitudinal stretching/compression on ice flow are ignored. The longitudinal deviatoric stresses are caused by extensional and compressional flow, which are likely to be important when describing the ice flow of valley glaciers (Egholm et al., 2011). This limitation is acceptable when modelling certain parts of ice sheets (althought it does not work well for ice streams or ice shelves), but can be more problematic for modelling mountain glaciers. The rugged terrain of mountain ranges presents some challenges for application of the SIA to alpine glaciers (Le Meur et al., 2004; Egholm et al., 2011). For example, this approximation neglects stress terms that are likely to be important for flow of mountain glaciers and ice sheets near their margins. Longitudinal stresses can be significant for glaciers in narrow or steep valleys. In an attempt to make the SIA more appropriate for modelling mountain glaciers, several studies have tried to account for longitudinal stresses. These adjustments either include parameterising the longitudinal deviatoric stresses (Hubbard, 2000) or solving for them directly in second-order SIA models (Baral et al., 2001; Egholm et al., 2011). Despite the fact that some of the assumptions that underlie the zeroth-order SIA are likely to be incorrect on alpine glaciers, comparisons between the SIA and full-Stokes models have shown similar results (Le Meur et al., 2004). For example, modelled ice extent and thickness are similar between the two approaches when modelling valley glaciers on gentle slopes (Leysinger Vieli and Gudmundsson, 2004). The major differences are in the small-scale velocity estimates generated by the two approaches. It has also been shown that the primary control on glacier evolution is usually the mass balance (Leysinger Vieli and Gudmundsson, 2004). For these reasons, a vertically-integrated 2-D ice-flow model based on the zeroth-order SIA can be used to model past extents of mountain glaciers (Plummer and Phillips, 2003; Kessler et al., 2006; Laabs et al., 2006). 44 Ice-flow models can be applied in 1, 2 or 3 dimensions. Flowline models (1-D) capture the basic ice dynamics along a centre flowline of a glacier. Flowline models require predefinition of a central flowline as well as parameterisation of valley geometry, where a valley cross profile is typically described using a trapezium or parabola. This also means that all tributary glaciers must be identified and parameterised as well. These parameterisations can lead to uncertainties in the results because, for example, the glacier width might evolve as ice thickens in a way that is not easy to parameterise. 1-D flowline models would also not be able to account for complex snow distribution caused by avalanching or snow drift, and are intrinsically less useful for comparing against geomorphological evidence, which is best presented in map-plan (2-D) form. I chose to use the 2-D, zeroth-order, vertically-integrated SIA model for several reasons; First, the valleys chosen in this thesis have relatively low angle slopes (5-15◦ ), for which the longitudinal stresses should be relatively low (Le Meur et al., 2004). Second, the resulting simulated glacier extents are in two dimensions (plan view) and are directly comparable with mapped moraines. Third, the vertically-integrated SIA is computationally fast, allowing for long integrations at relatively high resolutions. 3.4.1 Ice-flow model The modelled ice velocity is comprised of two components; sliding and deformation. Both sliding and deformation velocity in this work follow the equations suggested by Kessler et al. (2006). In this 2-D model, the vertical velocity profile is integrated, which ~ d ) is decreases computational time. The ice velocity caused by internal deformation (U calculated using: ~ d = 2 AH~τbn U 5 (3.23) where A is the coefficient of Glen’s flow law, set to 3.17e−25 Pa−3 s−1 , H is ice thickness (m), ~τb is the gravitational driving stress (~τb =ρice gH∇z), and n is Glen’s flow law exponent, set to 3. The sliding velocity follows the empirical formulation of Kessler et al. 45 (2006): τc ~ s = Uc e1− ~τb U (3.24) where Uc is a typical sliding velocity (m a−1 , Table 3.3) and ~τc is the characteristic gravitational driving stress (105 Pa, Kessler et al. (2006)). In Chapter 6 where the Cameron Glacier is modelled to its Holocene positions, Uc is made to vary with elevation to improve the fit in ice thickness to the height of late-Holocene lateral moraines. The Uc is set to 10 m a−1 for elevations greater than 1700 m asl in the model, and then increases linearly to 60 m a−1 for elevations lower than 1500 m asl. This imposed variability seems sensible because, at elevations below 1500 m asl, the valley floor flattens out and contains a layer of debris, which likely results in an increase in basal motion. In Chapter 4 where the Irishman glacier 13 ka position is simulated, the characteristic sliding velocity Uc is set to 20 m a−1 , following Kessler et al. (2006). Table 3.3 Ice-flow model parameter options and the range of published values for Glen’s flow law coefficient. Parameter Typical sliding velocity Symbol Uc Characteristic gravitational driving stress Ice density Glen’s flow law coefficient ~τc ρice A Other published values for Glen’s A Value 20 m a−1 10 - 60 m a−1 105 Pa 917 kg m−3 1e−17 Pa−3 a−1 3.17e−25 Pa−3 s−1 2e−24 Pa−3 s−1 1.3e−24 Pa−3 s−1 6.8e−24 Pa−3 s−1 4e−24 Pa−3 s−1 Source (Kessler et al., 2006) Chapter 6 (Kessler et al., 2006) (Paterson, 1994) This thesis (displayed in text) This thesis (different units) (Le Meur and Vincent, 2003) (Le Meur et al., 2004) (Kessler et al., 2006) (Cuffey and Paterson, 2010) The continuity equation (Equation 4.4) was used to evolve the glacier geometry through ~d + U ~ s ) were calculated at points offset time. Using a staggered grid, the ice velocities (U from where the ice thickness change was calculated. The flux gradients were then used to update the ice thickness using a forward explicit time-step. The evolving time-step was calculated using the stability criterion following Hindmarsh and Le Meur (2001). To minimise the error related to ice thickness at the boundaries, I applied an ice-flux correction. The correction was applied at every timestep to find areas where the ice flux 46 exceeded the ice thickness at a cell (Plummer and Phillips, 2003). This can occur at modelled glacier margins where the ice-free ground elevation of a cell is greater than the ice-surface elevation in an adjacent cell. In the correction, the ‘phantom ice’ mass was identified, removed from the glacier, and moved back to its original cell (Plummer and Phillips, 2003). 3.5 Model applications One overarching goal of this thesis is to understand the behaviour and sensitivity of four glaciers that deposited well studied moraine sequences in New Zealand. The combination of models that I used to simulate palaeo-glacier extents was designed for addressing the research questions (outlined in Section 2.7). The present-day topography, glaciology, and climatology are useful benchmarks for evaluating changes in past climate from today, which is why I used present-day topography and climatology in these simulations. In the steady-state model runs, a glacier length is specified by targetting a terminal moraine position and the model output provides an envelope of temperature and precipitation change combinations. The envelope can be further restricted by limiting the precipitation change to within a particular range (e.g. ±20 or 50%). Due to differences in glacier response time, the duration of each steady-state run was varied between 240 (Irishman glacier) to 300 years (Cameron Glacier). The steady-state simulations allow me to estimate the magnitude of past temperature and precipitation changes that are consistent with the Lateglacial and Holocene ice advances, thus addressing the first research question (i.e. quantifying past climate Section 2.7). Two of the studies in this thesis include transient simulations, where climate is varied and glacier extent continually adjusts during the model run. In Chapter 4 where the Irishman glacier was simulated, I forced the model with a previously published chironomid-based temperature reconstruction from Boundary Stream Tarn (Vandergoes et al., 2008), with no change in precipitation, to assess the fit between the modelled ice extent and the glacial geomorphic record. The Antarctic Cold Reversal (ACR) cold peak in the chironomid record caused an advance in the modelled glacier terminus to within 100 m of the well47 dated ACR moraine in Irishman Stream. In Chapter 6, I forced the glacier model with present-day interannual natural climate variability using ‘white noise’ (Oerlemans, 2000; Roe, 2011). This approach allowed me to examine the question of whether the Holocene glacier fluctuations in the Cameron valley were caused by climate shifts, or alternatively, from interannual stochastic variability in an steady climate. Table 3.4 Model information summary for comparison between the different studies. Location Time period Timestep DEM resolution DEM size MTD EBM Albedo mode Debris cover Ice-flow Ch. 4 Irishman Lateglacial monthly 25 m 241x201 Yes ELA Yes Ch. 5 Cameron 2009-2010 daily 25 m 161x169 Yes Yes time included - Ch. 6 Cameron Holocene monthly 50 m 180x184 Yes Yes ELA Yes App. A Whale Lateglacial monthly 100 m 71x101 Yes ELA Yes Each model includes various assumptions and uncertainties. Parameter values were varied in order to (1) provide minimal tuning to improve the fit between modelled ice extent and the glacial geomorphic record, (2) explore model sensitivity, and (3) estimate uncertainty in the palaeoclimate reconstructions. Parameter testing was also used in the snow mass transport and deposition model to achieve an optimal fit between measured and modelled snow distribution (Chapter 5) for the present-day. Sensitivity and parameter tests provide means of quantifying uncertainties in model results for palaeo-glacier simulations. In turn, these results allow a first-order assessment of uncertainty to be placed on the palaeoclimate estimates. 3.6 Summary Various models have been employed to study past glacier extents and the associated climate. I used three models in this thesis to simulate past glaciers: a distributed energybalance model, a snow transport and deposition model, and a 2-D ice-flow model. Each model contains a degree of approximation, and sensitivity tests allow the affect of these 48 uncertainties to be quantified in terms of the resulting palaeoclimate estimates. The intermediate level of complexity and fast computational execution of the models used here have allowed for hundreds of individual runs, which represent (albeit in a simplified way) relatively complex physical processes and a suite of different conditions. These models were chosen with the aims and research questions of this thesis in mind (Section 2.7). Together, the work within this thesis attempts to address some important questions about glaciers, moraines, and past climate. 49 50 Chapter 4 Evaluation of Lateglacial temperatures in the Southern Alps of New Zealand based on glacier modelling at Irishman Stream, Ben Ohau Range Abstract Climate proxy records from the middle to high latitude Southern Hemisphere indicate that a Lateglacial (15,000 - 11,500 years ago) climate reversal, approximately coeval with the Antarctic Cold Reversal (ACR), interrupted a warming trend during deglaciation. In New Zealand, some palaeoclimate proxy records indicate a cool episode during the ACR (ca 14,500 to 12,500 years ago), while others do not express a significant change in climate. Recently published moraine maps and ages present an opportunity to improve the palaeoclimate interpretation through numerical modelling of glaciers. We use a coupled energy-balance and ice-flow model to quantify palaeoclimate from past glacier extent constrained by mapped and dated moraines in the headwaters of Irishman Stream, a highelevation catchment in the Southern Alps. First, a suite of steady-state model runs is used to identify the temperature and precipitation forcing required to fit the modelled glacier to well-dated Lateglacial moraine crests. Second, time-dependent glacier simulations forced by a nearby proxy temperature record derived from chironomids are used to assess the fit with the glacial geomorphic record. Steady-state experiments using an optimal parameter set demonstrate that the conditions under which the 13000 year old moraine formed were 2.3 - 3.2◦ C colder than present with the range in temperature corresponding to a ±20% variance in precipitation relative to the present-day. This reconstructed climate change relative to the present-day corresponds to an equilibrium-line altitude of ca 2000±40 m above sea level (asl), which is ca 400 m lower than present. Time-dependent simulations of glacier length produce ice advance to within 100 m of the 13000 year old terminal moraine, indicating that the chironomid-based temperature forcing and moraine record 51 provide consistent information about past climate. Our results, together with other climate proxy reconstructions from pollen records and marine sediment cores, support the notion that temperatures during the ACR in New Zealand were ∼2 - 3◦ C cooler than today. 4.1 Introduction During the Lateglacial interval (∼15,000 - 11,500 years ago or 15 - 11.5 ka) in the middle to high latitudes of the southern hemisphere, a warming trend was interrupted by the Antarctic Cold Reversal (ACR, ∼14.5 - 12.5 ka (EPICA Community Members, 2006)). New Zealand is one of the few locations in the southern middle latitudes where Lateglacial moraines exist. Many examples of Lateglacial moraines in the Southern Alps of New Zealand have been identified (Porter, 1975b; Suggate, 1990; Fitzsimons, 1997), mapped (Birkeland, 1982; Barrell et al., 2011), and some have been dated (e.g., Turney et al., 2007; Kaplan et al., 2010; Putnam et al., 2010a). Transformation of glacier fluctuation records to climate changes is not a trivial process due to several non-climatic influences and is not readily achievable for some types of glacier (e.g., debris-covered or calving, (Winkler et al., 2010; Anderson and Mackintosh, 2012; Chinn et al., 2012)). Porter (1975a) estimated that the Lateglacial equilibrium-line altitude (ELA) in the Tasman River catchment was 500±50 m lower than modern, based on an accumulation-area ratio (AAR) analysis of reconstructed glaciers inside the Birch Hill moraine limit. Other proxy records rarely afford quantitative climate estimates. Besides glacial deposits, records that show a prominent cooling during the Lateglacial include sea surface temperature reconstructions from around New Zealand (Pahnke et al., 2003; Barrows et al., 2007a; McGlone et al., 2010), pollen records from upland locations in New Zealand (Burrows and Russell, 1990; McGlone, 1995; Newnham and Lowe, 2000; Turney et al., 2003; McGlone et al., 2004; Vandergoes et al., 2008), and chironomid-derived temperature reconstructions (Vandergoes et al., 2008). Less pronounced climate reversals within or overlapping with the ACR period have also been inferred based on other pollen (Newnham et al., 2007, 2012) and speleothem (Williams et al., 2005) records (see also Alloway et al. (2007)). In this paper, 52 we aim to improve our understanding of the Lateglacial climate event in New Zealand using glacier modelling techniques. 4.1.1 Glaciers as a climate proxy Glacier extent depends on the balance between accumulation and ablation (Oerlemans, 2001, 2005). Moraine positions in mountainous regions provide evidence of former glacier extents and may be used to infer past climates. Glacial geomorphic mapping of the central South Island of New Zealand (Barrell et al., 2011) has yielded a high-quality constraint on past ice extents, and several surface exposure dating studies provide ages for the culmination of Lateglacial ice advances (Ivy-Ochs et al., 1999; Kaplan et al., 2010; Putnam et al., 2010a). We present an application of a 2-D ice-flow and energy-balance model (EBM) to reconstruct temperature and precipitation values for the ACR. We apply these models to a glacier whose extent is delineated by moraines mapped by Kaplan et al. (2010) in the headwaters of Irishman Stream in the central Southern Alps of New Zealand and compare our findings to their AAR reconstructions and other climate proxy records. In experiment 1, we explore the different combinations of temperature and precipitation required to produce a steady-state glacier that most closely matches the moraine record, taking account of uncertainties associated with model parameter choices. In experiment 2, we use a chironomid-derived temperature reconstruction (Vandergoes et al., 2008) to drive our glacier model and test if this climate forcing produces ice extents compatible with the moraine record (Kaplan et al., 2010). 4.1.2 Study area The South Island of New Zealand is dominated by the high axial range of the Southern Alps, which impedes the prevailing westerly winds, resulting in a steep west to east precipitation gradient (Griffiths and McSaveney, 1983; Henderson and Thompson, 1999). Glaciers with fast response times to climate perturbations are abundant in the central 53 Southern Alps due to high snowfall rates and steep mountain slopes (Chinn, 1996; Fitzharris et al., 1999). We focus on moraines deposited in the headwaters of Irishman Stream (43◦ 59’30”S, 170◦ 03’00”E), located in the Ben Ohau Range (Figure 4.1) (McGregor, 1967; Birkeland, 1982; Kaplan et al., 2010). The dated glacial deposits lie at the edge of a cirque, referred to here as Irishman basin, that comprises the upper 2 km of Irishman Stream valley. The basin floor, for the most part, slopes gently (∼15◦ ) towards the southwest. Aside from the mapped and dated moraine sequence, there are several reasons why Irishman basin is well-suited for this study. First, Irishman basin is located <30 km to the southeast of the most glacierised section of the Southern Alps and is affected by the same regional atmospheric circulation, including the prevailing westerly winds (Kaplan et al., 2010; Chinn et al., 2012). Second, Irishman basin lacks complex terrain and has a shallow gradient, making our use of the shallow ice approximation equations an appropriate model choice (Leysinger Vieli and Gudmundsson, 2004). Third, the basin supported a single glacier with no tributaries, making it a very simple glacier system for ice reconstructions. Fourth, the present-day digital elevation model provides an appropriate topographic boundary condition for the numerical model because the basin is presently ice-free, and the visible bedrock knobs and outcrops on the valley floor indicate little sediment accumulation and hence minor topographic change since deglaciation. Moreover, the 15◦ down-valley slope of the basin floor precludes the former existence of large proglacial lakes and calving ice termini. Therefore, we assume glacier fluctuations in this valley directly relate to changes in atmospheric conditions and hence climate over this part of the South Island. There are several reasons why we believe that the Irishman basin glacier at 13 ka did not contain a significant surface debris cover. Although there are rock glaciers at the headwalls of Irishman basin and the basin has a veneer of rock debris, through which bedrock hills are visible, the Lateglacial moraine is not particularly large, especially the left lateral segment (Kaplan et al., 2010). These features suggest that debris cover of past glaciers in this valley was not sufficient to dramatically affect the albedo/energy balance. Also, present-day Glenmary Glacier, 12 km west of Irishman basin, is an example of 54 a small cirque glacier in a similar setting to that of the 13 ka Irishman glacier that is relatively free of debris-cover. There are no long-term climate data from Irishman basin. We characterised the present climate by interpolating between data from multiple low-elevation stations and measurements from a temporary climate station that we installed in Irishman basin (43◦ 59’30”S, 170◦ 2’57”E, 2010 m asl). The station operated over an 11 month period, March 2010 - January 2011, so February values are inferred to be the same as January temperature values to calculate an annual mean temperature of 1.2◦ C. Hourly temperatures from the temporary station were averaged into monthly mean temperatures (Table 4.1), with the warmest month at 6.2◦ C (January, 2011) and the coolest month at −6.1◦ C (June, 2010). Interpolated precipitation data suggest an annual mean of 1240 mm (Stuart, 2011) for this location. Although the annual total precipitation surface and monthly mean interpolated precipitation data are based on data from many long-term rain gauges located in the central Southern Alps, we use two nearby stations to provide an indication of precipitation variation at Irishman basin. Mean precipitation data from Twizel, 30 km to the south (Station numbers 4995, 4996, 4997), from 1973 to 1997 is ∼610 mm a−1 with a standard deviation of ±110 mm a−1 (18%) (National Institute for Water and Atmospheric Research (NIWA), Retrieved 2009-2011). Likewise, mean precipitation data from Aoraki/Mt. Cook village, 30 km to the north (Station numbers 18125, 4591, 4593), from 1930 to 2011 is ∼4100 mm a−1 with a standard deviation of ±900 mm a−1 (22%) (NIWA, Retrieved 2009-2011) (Twizel and Aoraki/Mt. Cook village locations are shown in Figure 4.1). We consider the characteristic present-day precipitation variability to be ∼20%, based on the standard deviation from these stations. 55 Figure 4.1 Location of the study area in central South Island, New Zealand. The Ben Ohau Range is outlined in green and Irishman basin is in the white box. Boundary Stream tarn (BST) is represented by the white dot southeast of Irishman basin and the white triangles are the locations of Aoraki/Mt. Cook village (north of the Ben Ohau Range) and Twizel (southeast of Ben Ohau Range). 4.2 Modelling glacier extent We used a simple iterative modelling approach to estimate palaeoclimate conditions based on modelled glacier extent. This strategy is similar to those employed in previous studies (Plummer and Phillips, 2003; Hubbard et al., 2005; Kessler et al., 2006; Laabs et al., 2006) except that our 2-D ice-flow model is coupled with a spatially-distributed EBM. Our EBM/2-D ice-flow combination is simple enough to simulate mapped ice extents quickly (steady-state runs take less than an hour on a desktop computer), yet is complex enough to account for shading, aspect, valley slope, and local climate. A digital elevation model was used to provide surface elevation throughout the model domain. This model was produced from 20 m interval topographic contours and spot heights on published 1:50,000-scale (NZMS260 series) topographic maps published by Land Information New Zealand. We simulated the 13 ka Irishman glacier at 25 x 25 m gridded horizontal resolution, with all simulations starting from an ice-free basin. The model domain covers a 5 x 6 km area centred over Irishman basin. 56 4.2.1 Input data Thirty years (1981 to 2010) of daily climate data were converted into monthly means for each grid cell in the model domain. We calculated climate data grids from several different sources depending on availability and reliability. First, daily climate data of relative humidity, solar radiation, temperature, and precipitation came from the NIWA CliFlo Database (NIWA, Retrieved 2009-2011). Relative humidity and solar radiation grids came from the virtual climate station network climate grid interpolations, available directly from the NIWA CliFlo website. The NIWA climate grids were not used for temperature because the data contain some biases in mountainous terrain (Anderson and Mackintosh, 2012). Instead, we created interpolated temperature surfaces using data from multiple, surrounding, low-elevation stations (NIWA, Retrieved 2009-2011). To make the interpolation surface, we took daily temperature data from each station (Tst in ◦ C), the monthly lapse rate ( dT in ◦ C km−1 ), and station elevation (zst in m) to create a ‘reference’ dz temperature (Tr in ◦ C) at sea level using: Tr = Tst − dT zst dz 1000 (4.1) Daily reference temperature surfaces were interpolated at sea level in the horizontal plane across the model domain (Tait and Zheng, 2007). Modelled temperature (T in ◦ C) in any grid cell was calculated using Equation (4.2), with the grid cell reference temperature, the month’s lapse rate, and the elevation value of the grid cell (z in m). Each modelled past temperature change (∆T ) is an additive change to the temperature (T ) in each grid cell and is uniformly applied across the domain: T = Tr + dT z + ∆T dz 1000 (4.2) We calculated and used variable monthly temperature lapse rates (Figure 4.2 and Table 4.1), which were determined by minimising the mismatch between Irishman temporary climate station temperature data and an interpolated temperature grid. This tempo57 rally varying lapse rate attempts to account for pervasive valley-scale temperature inversions which occur in winter in this region (Figure 4.2). Observed and calculated ‘clearsky’ solar radiation were used to calculate cloudiness, following Hock (2005) and Anderson et al. (2010). Second, reanalysis data from the National Centers for Environmental Prediction (NCEP) for present-day wind at the 850 hPa level (1981 - 2010) (Kalnay et al., 1996) were scaled to match observed wind speed and were applied uniformly across the model domain. Wind speeds are not modified for the complex topography and we assume wind speed distribution is uniform over the domain. These data were converted into monthly mean gridded wind speed. Further discussion of these datasets and the rationale for using them is provided in Anderson and Mackintosh (2012). The third source of input climate data came from an annual mean precipitation surface (Stuart, 2011) based on rain gauge data, from 1971 to 2000, spatially distributed around the Southern Alps to capture the steep precipitation gradient. The interpolated precipitation (Anderson et al., 2010) amount at each grid cell was calculated at monthly mean values from daily measured rainfall at lowland stations (NIWA). The interpolated data were guided by an annual mean precipitation surface (Stuart, 2011). Modelled past precipitation change (∆P ) is calculated as a percentage from present day, where present-day precipitation is ∆P =0% change, doubling present-day precipitation is ∆P =+100% change, and halving present-day precipitation is ∆P =−50% change. ∆P represents the change in precipitation reaching the site, which could be a result of regional changes in annual precipitation and/or changes in the amount of snow settling within Irishman basin as a result of increasing or decreasing snow transport. 4.2.2 The energy-balance model We developed a spatially-distributed EBM that uses the specified topography, monthly interpolated meteorological input, ice thickness at each monthly time-step, and prescribed climate perturbations to calculate annual mass balance at each grid cell. The resulting 58 Table 4.1 Measured monthly mean temperatures (Tmeas ) and calculated temperature ) for each month based on Irishman basin climate station data (2010-2011, lapse rates ( dT dz 43◦ 59’30”S, 170◦ 02’57”E, 2010 m asl). Annual mean measured temperature is ±1.2◦ C (assuming a value of 6.2◦ C for February) and the annual mean lapse rate is −5.4◦ C km−1 . Month Tmeas [◦ C] dT ◦ −1 dz [ C km ] Month Tmeas [◦ C] dT ◦ −1 dz [ C km ] ∗ Jan 6.2 -6 Jul -6 -3.5 Feb – -6∗ Aug -3.6 -5.5 Mar 4.9 -5.5 Sep -1.5 -6.5 Apr 3.1 -5.5 Oct 0.8 -5.5 May -1.1 -5 Nov 4.9 -5.5 Jun -6.1 -4 Dec 6 -6 was set equal to the January value because of a data gap. The February dT dz 8 6 Temperature (°C) 4 2 0 -2 -4 -6 -8 Mar Measured Temp Modelled Temp (variable lapse rate) Modelled Temp (lapse rate -6°C km-¹) Apr May Jun Jul Aug Sep Oct Nov Dec Jan Feb Month of the Year (2010-2011) Figure 4.2 Modelled versus measured temperatures at the location of a climate station in Irishman basin. Modelled temperature values (daily interpolated temperature, T , smoothed with a thirty-day running mean) were interpolated from climate stations at <800 m asl over the 2010-2011 interval. Measurements (black line) are compared to modelled temperatures that were calculated using a monthly varying lapse rate (red line) and a constant lapse rate (blue dashed). Using a constant lapse rate achieves a poor match, especially during winter months (Jun-Sep). 59 mass balance grid was then used as a boundary condition in our ice-flow model (Section 4.2.3). Several other workers have employed 1- and 2-D models based on similar sets of equations (e.g., Oerlemans, 1992; Plummer and Phillips, 2003; Hock, 2005). We followed an EBM scheme described in Anderson et al. (2010) and Anderson and Mackintosh (2012), except ours used a monthly timestep and updates mass balance once every 5 or 20 model years (transient and steady-state respectively). The model calculated the energy available for melt (QM ) using the following equation: QM = I(1 − α) + L ↓ +L ↑ +QH + QE + QG + QR (4.3) where I is the incoming shortwave radiation, α is the surface albedo, L ↓ is the incoming longwave radiation, L ↑ is the outgoing longwave radiation, QH and QE are the sensible and latent heat fluxes, QG is the ground heat flux, and QR is the heat supplied by rain. All heat exchanges are in units of W m−2 . Positive values of radiative and turbulent fluxes indicate a gain of energy for melting of the snowpack, and negative values indicate a loss of energy. QG was set at 1 W m−2 (Neale and Fitzharris, 1997) and QR was calculated assuming that precipitation is the same temperature as the air (Oerlemans, 1992). Incoming shortwave radiation (I) was comprised of direct and diffuse components following Oerlemans (1992). Cloudiness (see above) was included in the insolation calculations following Hock (2005) and Anderson et al. (2010). Insolation values for the steady-state simulations were calculated for 13 ka (Berger and Loutre, 1991; Eisenman and Huybers, 2006). The albedo parameterisation (α) was based on a background albedo profile, which shows albedo increasing with altitude and snow thickness, both dependent on the modelled ELA, following Oerlemans (1992). Longwave and shortwave radiation distribution were calculated to include the view field of the surrounding topography from the cell and cloudiness following Plummer and Phillips (2003) and Anderson et al. (2010). Turbulent heat fluxes (QH and QE ), which can make up half or more of the energy available for melt in maritime environments (Anderson and Mackintosh, 2012), were calculated using the bulk method (Klok and Oerlemans, 2002; Oerlemans and Grisogono, 2002) following Oerlemans (1992). The turbulent heat calculations used different roughness lengths 60 for snow and ice surfaces and the Richardson stability criterion was applied for stable stratification conditions (Oerlemans, 1992; Hock, 2005; Anderson et al., 2010). Snow accumulation occurred in the model when the temperature at each grid cell was at or below a snow/rain temperature threshold. The snow temperature threshold (Tsnow ) was set to 1◦ C (Anderson et al., 2006). We did not take possible surface debris cover effects into account in these palaeo-simulations, although this cover was likely insignificant when Irishman glacier filled much of the basin at the ACR extent as suggested by the present geomorphology and previous glacier reconstructions (Kaplan et al., 2010). Mass balance of the simulated glacier was calculated by subtracting the mass of possible melt from the mass of snow accumulation at each grid cell. The mass balance was recalculated every 20 model years in steady-state runs (every 5 model years in transient runs) while the ice-flow model ran. Coupling the models allowed the parameters sensitive to topography or the presence of ice (i.e., albedo, turbulent heat fluxes, shading, longwave radiation) to readjust as the ice thickness and extent evolved. 4.2.3 The ice-flow model The ice-flow model employed in this study is a 2-D shallow ice approximation model, with an explicit time-step, similar to Plummer and Phillips (2003) and Kessler et al. (2006). The model used the basic mass continuity equation: dH = M − ∇ · ~q dt (4.4) where H is the ice thickness, t is time, ~q is ice flux, and M is the mass balance, which was calculated using the EBM described in Section 5.2.2. The flux divergence was calculated using a 2-D finite-difference scheme. The vertically-averaged ice velocity from internal deformation was calculated using the shallow ice approximation (Paterson, 1994): 61 ~ d = 2 AH~τ n U b 5 (4.5) where A is the coefficient of Glen’s flow law, set to 1e−17 Pa−3 a−1 , ~τb is the gravitational driving stress (~τb =ρgH∇z), and n is Glen’s flow law exponent, set to 3. The sliding velocity follows the empirical formulation of Kessler et al. (2006): τc ~ s = Uc e1− ~τb U (4.6) where Uc is a typical sliding velocity (20 m a−1 ) and τc is the gravitational driving stress ~ s. (105 Pa) that results in U The continuity equation (Equation 4.4) was used to evolve the glacier geometry through ~d + U ~ s ) on a square grid that is offset from the time by calculating the ice velocity (U points at which ice thickness is known. The flux gradients were then used to calculate the updated ice thickness using a forward explicit time-step. The evolving time-step was calculated from the stability condition from Hindmarsh and Le Meur (2001). One way to quantify the error associated with the ice-flow model is to find the proportion of integrated mass balance for a modelled glacier in equilibrium relative to the total accumulation rate. The integrated mass balance over the glacier surface, for the Irishman glacier while it is in equilibrium (where ∆T = -2.5◦ C and ∆P = 11%), is -17.4 mm a−1 . The integrated total annual accumulation for the glacier is 815 mm a−1 . The mass balance divided by the accumulation total gives an estimate of the proportion of mass generated by the flow model (2.1%). The proportion is low, indicating that the ice-flow model is not a significant source of error for the palaeoclimate estimates. The assumptions that underlie the shallow ice approximation are less valid in parts of the model domain where bed slope and glacier aspect ratio are high (Le Meur et al., 2004), however, Irishman glacier overall has a low aspect ratio and a low ice surface gradient. We consider the primary control on glacier evolution in Irishman basin to be the mass balance (Leysinger Vieli and Gudmundsson, 2004), and thus the shallow ice approximation is 62 appropriate for our investigation. A summary of the parameter values used in the climate data interpolation, energy balance, and ice-flow models is given in Table 4.2. Most values were taken from published literature and value ranges in the sensitivity tests are used to calculate their effect on results. Table 4.2 Model parameter values. These values define our optimal parameter set, meaning they come from previously published studies optimal for modelling and/or are specific to the site location. Our sensitivity tests use these values, changing one value at a time to understand each parameter’s influence on our climate reconstruction. Parameter Name Snow albedo - ELA related Snow/rain temperature threshold Temperature lapse rates Roughness parameter for ice Roughness parameter for snow Typical sliding velocity Glen’s flow law coefficient Wind Maximum snow thickness 4.2.4 Symbol αsnow Tsnow dT dz Zice Zsnow Uc A Value 0.72 1◦ C see Table 4.1 0.004 m 0.001 m 20 m a−1 1e−17 Pa−3 a−1 Reanalysis 5 m w.e. Source Oerlemans (1992) Anderson et al. (2006) This study Anderson and Mackintosh (2012) Brock et al. (2006) Kessler et al. (2006) Paterson (1994) Kalnay et al. (1996) Steady-state simulations In experiment 1, incremental changes in temperature relative to modern day (∆T ) from −1.5 to −3.5◦ C, and in precipitation relative to modern day (∆P ) from −50% to +200% were applied uniformly across the year and domain to force the EBM/2-D ice-flow model from initial ice-free conditions. We assumed that the glacier was in or near equilibrium with the climate when it deposited the continuous moraine, and thus steady-state simulations represent a possible climate for a particular glacier size. The resulting steady-state ice extent was then compared to the moraine position, where a good match shows the distance between them being less than two grid cells (<50 m). These iterations were used to derive a ∆P -∆T curve, where combinations of ∆P -∆T on that curve produced a glacier extent that fitted the moraine position. Each steady-state simulation ran for 240 model years, which is greater than the maximum time for any ∆P -∆T combination for the 13 ka glacier to adjust to equilibrium starting from ice-free conditions. The EBM ran once every 20 model years to recalculate the glacier mass balance. These simulations used values from our optimal parameter set (Table 4.2). 63 In order to assess the impact of parameter choices on our results, the following parameters were systematically explored: Tsnow , αsnow , Zice , Zsnow , Uc , and A. In each case the model ran using the optimal parameter set with the exception of the parameter we were testing. These tests were done with a similar model setup to experiment 1, each running for 240 model years, changing ∆P -∆T combinations until the glacier reached within 50 m of the moraine. In the sensitivity tests, we created a new ∆P -∆T curve for each tested parameter (Figure 4.4). 4.2.5 Transient runs In experiment 2, transient model runs were forced using chironomid-derived mean summer temperature reconstructions from sediments at Boundary Stream tarn (BST, 44◦ 02’S, 170◦ 07’E, 830 m asl) located 7 km southeast of Irishman basin, that span the Lateglacial time (Vandergoes et al., 2008) (Figure 4.1). The age model at the Boundary Stream tarn (BST) site is derived from Bayesian modelling of 17 radiocarbon dates, calibrated using IntCal04 (Vandergoes et al., 2008). The chironomid record shows a period of temperature instability between 14.2 and 13.2 ka, including a maximum cool spike between 2 and 3◦ C below present values (Vandergoes et al., 2008). Particular strengths of the BST record are that chironomids are sensitive to mean summer temperatures (Dieffenbacher-Krall et al., 2007), the site lies close to Irishman basin, and the chironomid record is compatible with pollen-based palaeoclimate proxies from the same core. Although the chironomid temperature reconstruction describes mean summer temperature values, we apply it to our model as mean annual temperature changes. The sample specific error on the chironomidderived temperatures is ±1.4◦ C. Vandergoes et al. (2008) reconstructed summer temperatures from chironomids using two different modern-day chironomid-temperature calibration models and published both the raw and smoothed data from these models. They drew their conclusions (as much as ∼2 - 3◦ C cooling during the ACR) from the smoothed data of the partial least squares (partial least squares (PLS)-smooth) model. We used all four temperature reconstructions (PLS-raw, PLS-smooth, weighted average partial least squares (weighted average partial least squares (WAPLS))-raw, and WAPLS-smooth) to drive our EBM/2-D ice-flow model ∆T values. Each simulation was compared to the 64 moraine sequence dated by Kaplan et al. (2010) to determine the best fit out of the four temperature reconstructions. Temperatures were interpolated between the data points in each reconstruction from Vandergoes et al. (2008) using nearest neighbour interpolation. Precipitation was set to 0% change (present-day) for all transient runs. Each simulation ran for 4000 model years (15 - 11 ka) and the EBM ran once every 5 model years to recalculate the energy balance on the evolving ice surface including differences in insolation (I) due to changing orbital parameters during the period 15 - 11 ka. 4.3 Results 4.3.1 Experiment 1: Steady-state simulations of the ACR Irishman glacier extent The first experiment of our investigation involved carrying out a suite of equilibrium runs at set ∆P -∆T combinations to identify those that produced a good match with the Irishman basin 13 ka moraine position (Figure 4.3a+b). Figure 4.3b shows the fit between modelled ice thickness and the moraine record during a single simulation where ∆T is −2.5◦ C and ∆P is +11%. The modelled Irishman glacier at the 13 ka extent has an area of ∼1.8 km2 , volume of ∼0.15 km3 and a maximum thickness of 100 m. The modelled glacier response time to a 1◦ C decrease in ∆T is 102 model years, suggesting that the glacier would respond to centennial scale shifts in climate, rather than decadal. The relatively long response time is consistent with the low slope angle of Irishman basin and a relatively low mass turnover. Figure 4.3c shows the mass balance of Irishman glacier when ∆T is −2.5◦ C and ∆P is +11%. Mass balance ranges from 1 m w.e. accumulation at highest elevations to −2 m w.e. ice melt at the terminus. The modelled ELA (mean elevation where mass balance equals zero) in this simulation is 2000±40 m asl. The present-day ELA is above the elevation of Irishman basin and has been extrapolated to >2300 m asl from present-day 65 glacier snowlines (Porter, 1975b; Chinn, 1995; Lamont et al., 1999; Mathieu et al., 2009). To assess this estimate and to compare our palaeo-ELA with that of today, we applied the same climate data interpolation scheme and energy balance model to an idealised, cone-shaped and ice-covered mountain rising to 2500 m asl at the same location as Irishman basin. The modelled present-day ELA for Irishman basin derived using this method is ∼2400±40 m asl, only ∼100 m above the highest parts of the headwall of Irishman basin. Therefore, the ELA lowering that corresponds to the ACR event is ∼400±40 m. Figure 4.4 displays the ∆P -∆T curves using the optimal parameter set (Table 4.2) and variations in parameters. A similar glacier geometry to that shown in Figure 4.3b will form with any ∆P -∆T combination on the black curve in Figure 4.4, although the glaciers resulting from a larger precipitation increase are slightly thicker. With no change in temperature, a precipitation increase of 241% is necessary to simulate the moraine position (not shown in Figure 4.4). Varying the Tsnow and αsnow individually has the largest impact of all parameters on our climate results. Changes of ±1◦ C in the Tsnow alone influences the temperature reconstruction by ±0.3◦ C. Likewise, changing the αsnow alone by ±0.05 modifies ∆T by ±0.6◦ C. Despite the effects of parameter choice on our climate reconstructions, the overall conclusions do not change. Variations in other parameters have smaller influences on the ∆P -∆T curve. For example, varying the local summer insolation (I) between low (18 ka) and high (10 ka) values makes no direct difference in the modelled glacier extent or climate reconstruction. Table 4.3 lists the parameters, the values used in the sensitivity tests, and the ∆T needed to model Irishman glacier without a change in precipitation. The model-derived ∆T ranges from −3.3 to −2.1 with the optimal parameter set (Table 4.2) curve crossing at −2.7◦ C when ∆P is set to zero. We consider the temperature uncertainty (0.6◦ C) as half the range of temperature values (−3.3 to −2.1◦ C) shown in Figure 4.4. A smaller temperature change requires an increase in precipitation. For example, a smaller ∆T of −1.5◦ C requires a coincident increase in ∆P of 80% to achieve the 13 ka glacier extent. In summary, the model simulations in experiment 1 suggest a temperature reduction of 66 2.7±0.6◦ C sustained for a century is needed to produce the mapped extent without a change in precipitation. The optimal parameter set curve shows a ∆T of −3.2◦ C when ∆P is −20% and −2.3◦ C when ∆P is +20%. Our palaeoclimate estimates are most affected by changes in the Tsnow and αsnow parameters. The modelled glacier geometry and ELA are very similar to previous reconstructions (Figure 4.3d) based on geomorphology and AAR methods (Kaplan et al., 2010). Table 4.3 Parameter name, symbol, value, and the ∆T needed to grow the glacier to the 13 ka moraine without a change in precipitation (x-intercept), with the optimal parameter set yielding a ∆T of −2.7◦ C. Parameter Snow albedo Symbol αsnow Snow/rain temperature threshold (◦ C) Tsnow Temperature lapse rate (◦ C km−1 ) Roughness parameter for ice (m) dT dz Zice Roughness parameter for snow (m) Zsnow Incoming shortwave radiation (ka) I Typical sliding velocity (m a−1 ) Uc Glen’s flow law coefficient (Pa−3 a−1 ) A 4.3.2 Value 0.67 0.77 0 2 -6 0.002 0.008 0.0005 0.002 10 18 10 30 1e−16 1e−18 ∆T (◦ C) -3.3 -2.1 -3.0 -2.5 -3.0 -2.7 -2.8 -2.7 -2.8 -2.7 -2.7 -2.7 -2.7 -2.7 -2.7 Experiment 2: Transient run driven by chironomid-derived temperature reconstructions The second experiment of our study was to run four time-dependent simulations of glacier extent forced by the different BST chironomid-derived temperature reconstructions (Vandergoes et al., 2008). The resulting glacier length changes were then compared with the mapped and dated moraine sequence in Figure 4.3a. Figure 4.5 shows the PLS-smooth BST ∆T (present-day temperature is ∆T = 0◦ C) from 15 to 11 ka, and the lower plot shows the corresponding change in glacier length. The 13 ka moraine position (in terms of glacier length) is represented on Figure 4.5 by the orange diamond