Proglacial Lakes Control Glacier Geometry and Behavior During Recession

Ice‐contact proglacial lakes are generally absent from numerical model simulations of glacier evolution, and their effects on ice dynamics and on rates of deglaciation remain poorly quantified. Using the BISICLES ice flow model, we analyzed the effects of an ice‐contact lake on the Pukaki Glacier, New Zealand, during recession from the Last Glacial Maximum. The ice‐contact lake produced a maximum effect on grounding line recession >4 times further and on ice velocities up to 8 times faster, compared to simulations of a land‐terminating glacier forced by the same climate. The lake contributed up to 82% of cumulative grounding line recession and 87% of ice velocity during the first 300 years of the simulations, but those values decreased to just 6% and 37%, respectively, after 5,000 years. Numerical models that ignore lake interactions will, therefore, misrepresent the rate of recession especially during the transition of a land‐terminating to a lake‐terminating environment.


Introduction
Sustained glacier recession is driving the development of proglacial lakes at the termini of glaciers as meltwater accumulates within bedrock basins, behind moraine ridges and outwash fan heads, or is impounded by dead ice (Carrivick & Tweed, 2013). Glaciers in contact with a lake are receding more rapidly than their land-terminating counterparts, for example, in the Himalaya (Maurer et al., 2019;Tsutaki et al., 2019), in Alaska (Larsen et al., 2007;Willis et al., 2012), and in New Zealand (Chinn et al., 2012). Glaciers terminating in proglacial lakes can lose mass by several mechanisms in addition to melt from energy exchanges at the ice surface, namely, calving and subaqueous melting, collectively known as frontal ablation (Maurer et al., 2016;Sakai et al., 2009;Truffer & Motyka, 2016;Watson et al., 2020). Glaciers hosting proglacial lakes can, therefore, be partially decoupled from climatic forcing due to the effects of frontal ablation (King et al., 2018;Kirkbride, 1993), causing nonlinear responses in ice geometry and behavior (Benn et al., 2012). individual sites, usually spanning a single ablation season (e.g., Purdie & Fitzharris, 1999), and remote sensing investigations are limited to the last decade or so at most. Nonetheless, in order to have confidence in projections of glacier and ice sheet mass loss, it is critical to quantify changes in glacier and lake extent and behavior over longer timescales. It is imperative to understand if the trajectories of recent rates of mass loss from laketerminating glaciers (e.g., Dykes et al., 2011;King et al., 2019;Purdie et al., 2016;Quincey & Glasser, 2009) can be extrapolated into the future. A solution to this dearth of process data, and to the restricted timeframe from which it comes, is to assess the impacts of ice-marginal lakes on the evolution of glacial systems as represented in the Quaternary geological record.
The Quaternary geological record contains abundant evidence of the impacts of proglacial lakes on continental ice sheet geometry and behavior, and hence on climate and sea level (Krinner et al., 2004;Li et al., 2012;Stokes & Clark, 2004;Törnqvist et al., 2004). For example, studies of the Laurentide Ice Sheet during the Last Glacial Maximum (LGM) have shown that deglaciation was characterized by complex reorganizations in flow patterns and dynamics which, in part, have been attributed to proglacial lakes triggering and sustaining rapid ice flow (Utting & Atkinson, 2019). Many reconstructions of continental ice sheets during the last glacial cycle using ice flow models (e.g., Gregoire et al., 2016;Stokes et al., 2012) do not include proglacial lakes. A notable exception is the recent work of Hinck et al. (2019) who present an algorithm that efficiently identifies and fills isolated lake basins with the corresponding maximum water level and has been applied to ice sheet reconstructions of North America and the present day. However, with estimated maximum areal extents of several hundred thousand kilometers (e.g., 710,000 km 2 ; Glacial Lake Agassiz-Ojibway), ice sheet-lake reconstructions on this continental ice sheet scale are not directly applicable to mountain glaciers in regions such as the Himalaya, Patagonia, or New Zealand.
To date, no studies have isolated the effects of proglacial lakes specifically on glacier dynamics. Glacier sensitivity to the presence of a proglacial lake at its terminus and to the dynamics of that lake remains entirely unknown from the Quaternary record at the sub-ice sheet scale. Consequently, investigations are needed at the glacier scale to understand more about the rate of recession under the controls of proglacial lake development. The aim of this study is, therefore, to quantify the impact of an ice-contact proglacial lake on the geometry (grounding line position, ice margin shape, and ice thickness) and ice flow of a mountain glacier during recession. To achieve this, we used a well-resolved landform record as an empirical foundation for simple scenario modeling experiments.

Study Area
During the LGM, glaciers extended out from the Southern Alps icefield in New Zealand to form local piedmont lobes (Barrell et al., 2011;Golledge et al., 2012;James et al., 2019). The largest of which extended into and across the basin now occupied by Lake Pukaki and is referred to as the Pukaki Glacier (e.g., McKinnon et al., 2012;. At its maximum LGM extent the Pukaki Glacier was~85 km long and covered 1,300 km 2 (Figure 1; Porter, 1975). Early geological and geomorphological work relating to the Pukaki Glacier is summarized by Barrell and Read (2014), but of most significance is that the Pukaki valley contains an exceptionally well-preserved geomorphological record of lateral and terminal moraines (Barrell et al., 2011). The moraines trend north to south and the ridges parallel the eastern shore of Lake Pukaki and have been dated in many geochronological campaigns (e.g., Doughty et al., 2015;Kelley et al., 2014;Schaefer et al., 2006Schaefer et al., , 2015Strand et al., 2019). Those moraines most relevant to this study are an inner line at the southern end of the lake, with a mean surface exposure dating (SED) age of 17.6 ± 1.2 ka, thus marking a deglacial age after the onset of the LGM termination (Darvill et al., 2016; Figure 1) and the Birch Hill moraines at the northern, up-valley end of the lake which have a mean SED age of 13.1 ± 0.5 ka (Darvill et al., 2016; Figure 1) representing either a still stand or a readvance of the ice margin during the Antarctic Cold Reversal (ACR; . It has been argued that the time gap of~5,000 years and the large valley distance of >25 km between these sets of moraine clusters is evidence for an abrupt glaciation termination with virtual ice collapse, in tandem with Antarctic warming (Schaefer et al., 2006).

Model Description and Domain
In this study, we use the vertically integrated ice flow model BISICLES (Cornford et al., 2013). BISICLES was chosen for its efficient and accurate representation of the dynamics of marine grounded ice sheets and for its process-representation permitting ice shelf formation, grounding line migration, and ice surface lowering. BISICLES has previously been applied to simulations of both contemporary (e.g., Antarctica, Cornford et al., 2015;Berger et al., 2016) and paleo ice sheets (e.g., The British-Irish Ice Sheet; Gandy et al., 2018Gandy et al., , 2019, as well as to smaller icefields (e.g., Patagonia). The dynamical equations in BISICLES fall into type "L1L2" of hybrid ice sheet modeling approaches (Hindmarsh, 2004), where the longitudinal stresses are treated as depth dependent, and are included in the computation of stresses driving the ice flow (Schoof & Hindmarsh, 2010). This approach includes elements from both the shallow ice approximation and shallow shelf approximation, which is more appropriate for modeling fast flowing and floating ice than either the shallow ice approximation or shallow shelf approximation alone (Bueler & Brown, 2009).

Experimental Design
We "spin-up" the model to produce a steady state ice geometry representative of the end of the LGM at 18 ka ( Figure S2 in the supporting information) as given from field evidence (Barrell et al., 2011) and by other modelled reconstructions (e.g., Golledge et al., 2012;James et al., 2019). Recession of two scenarios were then simulated: a land-terminating Pukaki Glacier ("LAND") and a lake-terminating Pukaki Glacier ("LAKE"). LAND was forced by an idealized climate only, whereas LAKE was forced by the same idealized climate in addition to frontal ablation (subaqueous melting and calving). Our scenario experiments were informed by geochronological data from the moraines at Pukaki but were not best fitted to them due to the lack of evidence of ice margin position in the basin between the LGM and the ACR. The position of the moraines separated by >25 km (Figure 1) enabled Porter (1975) to estimate an equilibrium line altitude (ELA) for the LGM and Birch Hill ice extents of 1,225 and 1,600 m above sea level, respectively. Assuming that the glacier was in steady state at both times, the difference in ELA and the mean exposure ages of the moraines, representing~5,000 years between formation, enabled us to calculate a linear rate of ELA rise. Climate forcing is therefore specified in our model via a linear increase in ELA of 0.05 m per year, or for clarity 50 m ELA rise per millennium. To capture deglaciation from the end of the LGM to up-valley of the Birch Hill limits, we ran both LAND and LAKE simulations from the spin-up extent ( Figure S2) for 10 kyr.

Boundary Conditions and Parameterization
Specifying an idealized climatic change, initiating the scenarios with an equilibrium ice extent initial condition and keeping all other boundary conditions and parameters identical within both sets of model simulations, enabled us to isolate the internal mechanisms of recession, specifically assessing the relative role of frontal ablation on ice dynamics. Information about the model initial conditions such as bed topography is discussed in Text S1 and Figure S1. Key model parameters are summarized in Table S1. We used a linear viscous sliding law (m = 1), which accounted for sliding that did not require basal melt. Our choice of calving and subaqueous melt parameter values in LAKE were informed by our sensitivity testing (Text S4 and Figure S3), which itself was based on values from contemporary observations reported in the literature (Table S2). We ran 10 simulations of LAKE: 5 prescribing a different subaqueous melt rate (0, −1, −10, −50, and −100 m a −1 ) and 5 prescribing a different calving rate (10, 50, 100, 500, and 1,000 m a −1 ;  Table S3) each time. We report LAKE as the calculated mean output value from these simulations. Our sensitivity testing revealed that subaqueous melt rate had a relatively larger effect on glacier geometry and behavior, but that calving rate did not. The subaqueous melt flux was masked only to floating ice. Calving into the lake was modeled according to a crevasse depth criterion (Benn et al., 2007;Nick et al., 2010). Surface crevasses opened to the depth where tensile stress and water pressure are balanced by ice overburden pressure (Weertman, 1973), and calving took place when the crevasses reached the waterline. Water pressure was set by assuming that crevasses filled with up to 100 m of water. Noting an absence of field data on fluctuating lake level, the lake level was held constant throughout the LAKE simulations.

Geophysical Research Letters
There are several differences to consider when applying an ice flow model, designed to simulate marine-terminating environments, to lacustrine settings (Text S3). To account for the differences between marine and lake-terminating glacier termini, we made the following specifications: (i) the density of freshwater (1,000 kg m −3 ; Table S1) was used to compute buoyancy as opposed to the density of saline water; (ii) smaller calving and subaqueous melt rates were used than are usually prescribed within ice sheet models (Tables S2  and S3); and (iii) water temperature was accounted for within the subaqueous melt flux component but held constant because thermal stratification is assumed unimportant on millennial timescales (Carrivick & Tweed, 2013).

Results
We computed grounding line recession of the Pukaki Glacier over tens of kilometers and examined the response of glacier geometry (grounding line position, ice margin shape, and ice thickness), ice flow, and rate of recession to the perturbation of terminus environment (LAND versus LAKE). Both simulations caused recession of the Pukaki Glacier terminus from the lip of its overdeepened basin northward past the Birch Hill moraines (Figure 1) within 6 kyrs ( Figures S4-S6), thus in broad agreement with the rate of recession interpreted from the glacial geomorphology (e.g., . However, in contrast to the geomorphological record, our simulations suggest the time-transgressive nature of ice margin recession. Substantial differences in grounding line position (terminus position in LAND), ice thickness, and velocity were simulated between LAND and LAKE from 16.5 to 14.5 ka (Figures 2 and 3). We specifically focus on comparing glacier geometry and behavior at the grounding line although the effects of LAKE on ice thickness and velocity were propagated up-glacier ( Figure 2).

Ice Margin Geometry
The Pukaki Glacier in LAND was grounded across its entirety throughout the duration of the simulation (Figures 2b and S4). In LAKE, the glacier was wedge shaped and floating where the grounding line remained 2 km behind the terminus (Figure 2b). The width of the glacier in LAND was much narrower in comparison to that within LAKE, which spanned the width of the overdeepened basin ( Figure S5) as an ice cliff. Throughout recession, the glacier in LAND became narrower ( Figure S5), whereas in LAKE it remained at the same width throughout recession. The maximum thickness of the Pukaki Glacier in both model simulations was 1,300 m (at 18 ka), occurring in the trunk of the glacier in the narrowest part of the valley ( Figure S2). Ice thickness in LAND was thinner, increasingly so toward the terminus (Figures 2 and 3), and varied between 15 and 74 m. During the first 100 years after the glacier receded into the basin in LAKE, progressively thicker ice floated; the thickness of ice at the grounding line in LAKE increased substantially up to 345 m at 16.9 ka (330 m thicker than LAND; Figure 3) and then steadily declined throughout the rest of the simulation. Lake effects contributed 96% of the total ice thickness after 100 years (at 16.9 ka) and diminished to 48% by the end of the simulations.

Ice Margin Position
During the first 1,000 years (from 18 to 17 ka), the terminus receded 2 km in both simulations (Figures 2b  and 3) but remained grounded on the lip of the overdeepened basin (Figure 2b). This recession was land-terminating and so driven by the rising ELA only. With continued recession, the Pukaki Glacier in the LAKE and LAND scenarios diverged from each other. The maximum divergence occurred from 17 to 16.9 ka at the time of entering the over-deepened basin. The position of the ice margin in the two scenarios then slowly converged again over the next 1,500 years (by 14.5 ka; Figure 3). The rate of recession slowed as the glacier receded through its over-deepened basin (Figures 2b and 3). The difference between simulations in position 10.1029/2020GL088865

Geophysical Research Letters
of the grounding line in LAKE and the terminus in LAND was greatest at 16.5 ka where the grounding line had receded 14 km by 16.5 ka in LAKE compared to 3 km of terminus recession during the same period in LAND (Figures 2b and 3). At 16.5 ka, the grounding line in LAKE had receded nearly 6 times further than the terminus of LAND. Lake effects contributed 76% of total grounding line recession after 100 years (at 16.9 ka). This increased to 82% of cumulative grounding line recession after 300 years before diminishing to 6% by the end of the simulations.

Ice Velocity
Maximum ice velocities of 141 m a −1 occurred in both simulations at 17.5 ka in the main trunk of the glacier where the ice was thickest (~48 km up-glacier where the valley narrows just north of the confluence of the Jollie valley along the X-Y transect in Figure 2a). Mean ice velocities at the terminus in both simulations prior to 16.5 ka were between 10 and 20 m a −1 when the ice margin was grounded on the lip of the overdeepened basin. The greatest differences in ice velocity between LAND and LAKE were at the terminus ( Figure 2). Grounding line velocities in LAKE were amplified above LAND terminus velocities by 1 order of magnitude (Figure 3). Low mean velocities varied in LAND between 11 and 47 m a −1 , whereas velocities were relatively increased at the calving ice margin of LAKE (Figure 2b). The lowest velocities in LAKE occurred at the lateral ice margin where the ice was grounded on the valley sides (Figure 2a). Figure 3  (Table S3) in LAKE simulations, and the blue line represents the mean of these output values.

Geophysical Research Letters
displays 700 years of acceleration at the grounding line in LAKE, showing steadily increasing grounding line velocities until reaching a peak of 109 m a −1 at 16.3 ka. The biggest difference in velocity between the LAND terminus and the grounding line of LAKE was 88 m a −1 , which occurred at 16.8 ka. Lake effects contributed 87% of the total ice velocity after 100 years (at 16.9 ka) and diminished to 37% by the end of the simulations.

Discussion
Our numerical modeling illustrates and explains the substantial effects that a lake at the terminus of a glacier can have on the geometry (grounding line position, ice margin shape, and ice thickness), ice flow, and rate of recession when compared to the same glacier terminating on land and under the same climate. These effects can be summarized as being more rapid ice margin recession, ice surface drawdown, and ice velocity acceleration (Figures 2b and 3), and in that regard they confirm prevailing conceptual models (e.g., Carrivick & Tweed, 2013). However, in this study we have quantified the spatiotemporal nature of these effects, and we are able to shed light on the driving mechanisms responsible for these differences.

Glacier Geometry and Behavior
Differences in glacier geometry and behavior between the simulations became more marked as the lake-terminating glacier encountered progressively deeper water within its overdeepened basin and frontal ablation became proportionally more effective (Figure 2). The divergence between the two simulations, that is, the effect of the proglacial lake, diminished once the glacier terminus was receding up a normal-grade bed slope (Figures 2b and 3). The initial rapid recession indicates that the lake-terminating glacier was abruptly forced out of equilibrium with the prescribed climate warming. Once the effect of the lake was no longer pronounced, the glacier entered a period of adjustment to re-equilibrate with the prescribed warming.
The land-terminating scenario displayed a linear rate of recession after 16.5 ka, whereas the lake-terminating scenario displayed an initial acceleration of recession by 16.9 ka with a near-constant rate of recession thereafter (Figure 3). After their initial differences, the rates of recession between the two scenarios were similar. We identify a "period of collapse" in the lake-terminating scenario from 17 to 16.5 ka. The timing of collapse was simulated to occur when the ice margin receded from the lip of its overdeepened basin, down a reverse-grade bed slope, and into deep water (>200 m; Figures 2b and S1). The early amplified recession during this period is explained by recession of the grounding line into its ice-contact proglacial lake that was both expanding laterally and deepening as the glacier receded northward into its basin (Figure 3). The reverse gradient of the bed slope persists for~4 km (Figures 2b and S1) of recession. This, and the transition of the terminus to a position of a shallowing and then of a normally graded bed slope (Figure 2b), explains why the effect of amplified recession does not increase through time.
We observe that the faster recession of the Pukaki Glacier before 16.5 ka in LAKE is related, at least in the model, to a phenomenon known in tidewater glaciers and marine ice sheets as marine ice sheet instability (Katz & Worster, 2010;Schoof, 2007;Thomas & Bentley, 1978;Weertman, 1974). All else being equal, recession of the grounding line along a reverse-grade slope (deepening toward the center of the glacier) results in faster flow and further recession. This effect can potentially be countered by buttressing (Gudmundsson, 2013) or spatially variable friction (Brondex et al., 2017). However, both are unlikely conditions in ice-marginal lakes. We propose that the processes and feedbacks involved in marine ice sheet instability extend beyond the marine realm to glaciers that terminate in water in large overdeepened basins. Indeed, Larsen et al. (2015) suggested that the presence of large proglacial lakes that occupy substantial overdeepenings may lead to an ice margin geometry that is more akin to tidewater glaciers; thus, similar feedbacks may operate to induce ice front recession and ice loss. The duration and magnitude of ice loss from lake-terminating glaciers will therefore depend on (i) the size of each glacial overdeepening, with glaciers only reaching equilibrium once the terminus can recede to an altitude above the lake they host, and (ii) ice thickness because that determines the threshold for flotation.

Wider Implications
The timing of recession of the Pukaki Glacier is thought to be consistent with the onset of temperature and atmospheric CO 2 increases indicated by Antarctic ice cores (Schaefer et al., 2006). However, we propose that this atmospheric warming likely only initiated glacier recession from its LGM position, which forced the terminus into its overdeepened basin and eventually into deep lake water. After which, our modeling suggests 10.1029/2020GL088865 Geophysical Research Letters that both frontal ablation processes and climatic warming contributed to local glacier collapse. Once the Pukaki Glacier became lake terminating, its terminus recession was likely partially decoupled from climate. Such hypotheses have been previously proposed in New Zealand for other glaciers (e.g., Shulmeister et al., 2019;Sutherland, Carrivick, Evans, et al., 2019; and have been confirmed for glaciers elsewhere by empirical data, for example, from Patagonia (Bendle et al., 2019). This partial decoupling of glacier behavior from climate is potentially important for interpreting the evolution of New Zealand mountain glaciers, particularly over decadal timescales.
The large changes observed from the fine temporal resolution (100 years) in our model (Figure 3) reveals variability between the LAND and LAKE scenarios on decadal to centennial scales but not over millennia. The rapid acceleration in recession followed by a near-constant rate of recession, and the behavior captured within our lake-terminating simulation has marked similarities with observations from contemporary glaciers. Rapid changes to lake-terminating glaciers have been shown to occur over two to three decades in the central Himalaya (Basnett et al., 2013;King et al., 2019;Song et al., 2017). King et al. (2018) specifically showed that Himalayan glaciers can undergo periods of ice acceleration and enhanced mass loss with the onset and continuation of lake formation, followed by ice deceleration and mass loss stabilization as the glacier recedes out of the overdeepening. Based on this empirical evidence and our numerical modeling, we suggest that this spatiotemporal pattern is common and virtually all lake-glacier interactions will evolve in this way, where the magnitude of a lake effect and the time duration of those effects are proportional to the size/depth and elongation of the lake basin, as well as to individual ice geometries and ice dynamics.
The persistence of proglacial lake effects upon ice dynamics requires more detailed study. Specifically, decadal-to centennial-scale projections will require ice sheet and mountain glacier modeling that incorporate proglacial lake evolution. With lake-terminating glaciers becoming more numerous and larger (e.g., Carrivick & Quincey, 2014), the likely increasing influence of frontal ablation processes must be incorporated into ice dynamics models with account for fluctuating lake levels driven by contributions to the meltwater flux. Future work should also consider the importance of several other mechanisms that could further enhance recession of a lake-terminating glacier in comparison to a land-terminating glacier, such as convective warming of the ice by thermal warming of lake water, and wind-driven wave erosion (Mallalieu et al., 2020).

Conclusions
This study has (i) assessed the timescale over which the impact of lake processes influence mountain glacier dynamics and (ii) quantified the importance of lake processes relative to climatic forcing. We have shown that a glacier terminus in contact with a lake receded a maximum of more than 4 times further and by up to 8 times faster than a land-terminating ice margin under the same climate warming. This acceleration was short-lived, however, relative to the overall timescale of our simulations, and these maximum lake effects occurred for <10% of the simulation when the lake initially made contact with the ice margin. After which, the two scenarios followed similar trajectories. Our results suggest that 82% of cumulative grounding line recession, 87% of total velocity, and 96% of total ice thickness can be attributed to lake effects. These are maximum values, and they occur during the early stages of recession (within a few hundred years) as the ice margin recedes down a retrograde slope, before diminishing in effect to 6% of grounding line recession, 37% of total velocity, and 48% of total ice thickness by the end of the simulations. We attribute the divergence in behavior of the lake-terminating scenario, relative to the land-terminating scenario, to the combination of a reverse-grade bed slope and frontal ablation resulting in increased ice flow rates at the grounding line and driving stress as the control. When meltwater is impounded within the overdeepened basin, the glacier can (i) float, thereby permitting basal melt, and (ii) calve, thereby permitting mechanical mass loss. Both mechanisms act independent of climate-driven surface melting and produce accelerated ice flow causing compression of the glacier tongue thus steepening the ice surface slope and increasing longitudinal driving stress. The glacier progressively thins and accelerates until the bed slope becomes normally graded, and the ice surface slope matches the bed slope. Although these processes have been conceptualized, our study has added spatiotemporal quantification.

Geophysical Research Letters
Climate warming is important in order to initiate glacier recession and the development of an ice-contact proglacial lake. However our numerical modeling highlights the transition from a land-terminating to a lake-terminating environment can dramatically perturb glacier geometry and behavior. This partial decoupling from climatic forcing is particularly important over decadal to centennial timescales. We contend that ice sheet outlets or mountain glaciers that develop ice-marginal lakes will experience accelerated mass loss, where the absolute spatiotemporal magnitude of that mass loss relates to the depth and elongation of the associated overdeepening, the lake water depth, and glacier characteristics such as ice thickness and velocity.

Data Availability Statement
The digital elevation model used is available from Land Information New Zealand (2019)