Toward Numerical Modeling of Interactions Between Ice-Marginal Proglacial Lakes and Glaciers

Global climate change is evidently manifest in disappearing mountain glaciers and receding and thinning ice sheet margins. Concern about contemporaneous proglacial lake development has spurred an emerging area of research seeking to quantitatively understand lake - glacier interactions. This perspectives article identifies spatio-temporal disparity between the coverage of field data, remote sensing observations and numerical modeling efforts. Throughout, an overview of the physical effects of lakes on glaciers and on ice sheet margins is provided, drawing evidence together from very recent and high-impact studies of both modern glaciology and of the Quaternary record. We identify and discuss six challenges for numerical modeling of lake - glacier interactions, namely that there are meltwater exchanges between glaciers and ice-marginal lakes, lake bathymetry and glacier bed topography are often unknown, lake - glacier interactions affect the longitudinal stress regime of glaciers, lake water temperature affects glacier melting but is very poorly constrained, the interactions persist with considerable spatio-temporal variability and with boundary migration, and data for model parameterization and validation is extremely scarce. Overall, we contend that numerical modeling is a key frontier in the cryospheric sciences to deliver process understanding of lake - glacier interactions.

In recent Quaternary studies, ice-marginal lakes that formed immediately following the Last Glacial Maximum (LGM) have been investigated in great detail; across the Northern Hemisphere in Europe and Asia (e.g., Batchelor et al., 2019;Emery et al., 2019;Utting and Atkinson, 2019;Dalton et al., 2020;Herget et al., 2020;Winsemann and Lang, 2020;Turzewski et al., 2020), in Greenland Adamson et al., 2019), in South America (e.g., Hall et al., 2019;Martin et al., 2019;Thorndycraft et al., 2019;Davies et al., 2020;Mendelová et al., 2020), and the Southern Alps of New Zealand (Sutherland et al., 2019a, b). Some of these lakes have been implicated in destabilizing ice sheets (Colman, 2002) and are thought to have been a control on ice stream onset and dynamics (Stokes and Clark, 2004;Perkins and Brennand, 2015). They have been associated with glacier surges (Clark et al., 1996;Ponce et al., 2019) and with meltwater pulses into the North Atlantic (Clark et al., 1996;Bauer et al., 2004). However, with a few exceptions (e.g., Clarke et al., 2009) these implied effects remain largely untested within numerical models.
Ice-marginal lake development must be understood if mountain glacier and/or ice sheet evolution; dynamics, mass loss, and meltwater contributions and sea level rise are to be accurately quantified. Unfortunately, field data are limited in spatial and temporal coverage. Unraveling cause and effect in remote sensing studies is difficult; is meltwater from diminishing glaciers making the lakes bigger, or are lakes accelerating glacier shrinkage, or is it (most likely) a combination of both and what is the balance between these sets of processes? Furthermore, the vast majority of glacier and ice sheet models to date, whether applied to modern or to Quaternary systems, either completely neglect ice-marginal lake effects, or else have an exceptionally trivial representation of them. This perspective article overviews the physical effects of ice-marginal lakes on glaciers and then discusses the key concerns for future development of numerical models.

PHYSICAL EFFECTS OF LAKES ON GLACIERS AND ICE SHEET MARGINS
Ice-marginal lakes interact with glaciers via a series of thermomechanical processes Carrivick and Tweed, 2013;Tweed and Carrivick, 2015;Truffer and Motyka, 2016). These processes act in addition to climate-driven ablation effects, but can interact with those via several feedbacks. Temporary destabilization of an ice-margin, partially decoupling a glacier response from climate forcing (Kirkbride, 1993), can occur via calving, causing draw-down of ice as ice surface slopes are steepened ( Figure 1). Furthermore, as a hydraulic connection is established between a subglacial drainage system and lake water it thereby reduces bed friction and reduces longitudinal stresses. Increased glacier flow speeds can be expected, with a reduction in back stress and decreased bed friction. The result is dynamic thinning (e.g., Trüssel et al., 2013;Tsutaki et al., 2019) introducing a positive feedback of decreasing effective pressure, enhanced basal ice motion and further enhanced flow velocity .
Mass loss at lake-terminating margins of glaciers and ice sheets can be attributed to the effects of "frontal ablation," which is a term used to include calving and subaqueous and subaerial frontal melt (Cogley et al., 2011). In contrast to marine glacier margins, there have been few attempts to quantify frontal ablation at lake-terminating ice-margins. Several studies have inferred subaqueous melt to be a negligible component of ablation (e.g., Boyce et al., 2007;Trüssel et al., 2013), whereas others have considered it to be a major component of mass loss and a rate-controlling process on calving, particularly where lake temperatures are relatively warm (e.g., Haresign and Warren, 2005;Sugiyama et al., 2016;Minowa et al., 2017) and/or the terminus is slow-flowing with limited crevassing (Kirkbride and Warren, 1997;Röhl, 2006). Subaqueous melt rates are a key control on the propagation of thermo-erosional notches at the waterline, and thus melt-undercut driven calving losses (Diolaiuti et al., 2006;Röhl, 2006). A recent analysis of glacier calving into an ice-marginal lake found associations between calving style and magnitude with presence/absence of lake ice cover and lake water depth (Mallalieu et al., 2020). Where lakes freeze or have a floating terminus, a buttressing effect (Geirsdottir et al., 2008;Tsutaki et al., 2013;Mallalieu et al., 2020) can be imparted on the glacier, rather like that of an ice mélange in a marine setting.

EFFECTS OF SUDDEN LAKE DRAINAGES
Ice-marginal lakes can suddenly drain due to failure of a moraine or ice dam. Typical triggers for such failure include mechanical breakage of unconsolidated moraine material or lake water pressure (due to water depth) exceeding ice overburden pressure (Tweed and Russell, 1999). Studies on modern glaciers have shown that when an ice-marginal lake drains, a local portion of the adjacent glacier is near-instantaneously de-buttressed promoting calving (e.g., Sugden, 1985;Mallalieu et al., 2020), localized crevassing, and an abrupt increase in horizontal velocity (e.g., Walder et al., 1996;Furuya and Wahr, 2005;Tsutaki et al., 2013). These effects have been detected with in situ measurements using global positioning systems (GPS) (e.g., Walder et al., 2006;Tsutaki et al., 2013), satellite interferometry (e.g., Furuya and Wahr, 2005), or via processing of multiple oblique field photographs (e.g., Mallalieu et al., 2017Mallalieu et al., , 2020. Vertical displacements of an ice surface have been interpreted to signify passage of energetic and pressurized lake water through/ beneath a glacier (e.g., Walder et al., 1996;Anderson et al., 2005;Tsutaki et al., 2013). Lake depth and longitudinal expansion of the lake are key to explaining glacier response to sudden lake drainage events, but are themselves a product of glacier dynamics, thereby complicating process understanding (Kienholz et al., 2020). Many ice-dammed lakes drain multiple times establishing a filling-draining cycle (Clague and Evans, 2000;Carrivick and Tweed, 2019;Godbout et al., 2019), which further complicates attributions of cause and effect. During lake (re)filling, abrupt changes in glacier surface elevations have been associated with the formation of faults and fractures forced by englacial water, as evidenced by fluctuations in shallow borehole water pressures, and changes in radar internal reflection power (Bigelow et al., 2020).
Whether outburst floods from glacial lakes need to be included in ice sheet evolution models remains unconsidered, despite recognition in the geological record of their widespread influence on subglacial hydrology, especially in association with ice streams, for example, around the periphery of Antarctica (Kirkham et al., 2019;Larter et al., 2019) and in northern Greenland (Keisling et al., 2020).

DISCUSSION OF CONCERNS FOR NUMERICAL MODELING
The spatio-temporal disparity between modern process-based field data and knowledge gained from the geological record mean that numerical modeling is necessary to deliver processunderstanding, specifically on the relative importance of contributions of ice-marginal lakes to overall glacier ablation and dynamics. The space and time scales over which such nonlinear feedback processes might be important need unraveling. A quantitative consideration of ice-marginal lake effects should be considered as a key frontier in glaciology. Six major issues are identified from the literature and are discussed in the next sections.

Meltwater Exchanges
As glaciers lose mass, retreat and thin, they release meltwater that can become ponded behind moraine ridges and within other topographic depressions. This is especially the case when a glacier rests on a retrograde slope (Nick et al., 2009), such as a glacial overdeepening, or on the periphery of an ice sheet covering an isostatically-depressed region, such as west Greenland (Morlighem et al., 2017). Until recently, representation of feedback between glacier meltwater generation and ice-marginal lake evolution remained beyond the capability and capacity of mountain glacier or ice sheet models.
Depressions in topography (see Figure 2) act as a sink for many algorithms of water flow routing (Zhu et al., 2006). Socalled "flood-fill" algorithms can be used to determine the extent and depth of lakes from topography data. However, computational time can become a limiting factor when applying flood-fill algorithms to large geographical extents and long timescales (for example, in the case of palaeo-ice sheets). Berends and van de Wal's (2016) model can be completely coupled to glacio-isostatic adjustment (GIA) models and thus the mass of lake water is added to the mass of ice for GIA/sea level calculations, accounting for bedrock deformation and geoid perturbation. These deformation effects both influence the volume of the lake by deepening the basin and through raising the water level. In contrast, Hinck et al. (2020) identify lake basins and the corresponding maximum water level for a given ice sheet and topography reconstruction without including advanced flow routing techniques. They applied their model to reconstructions of North America between the LGM and the present day, successfully reconstructing Lake Agassiz, Lake McConnell and the predecessors of the Great Lakes. Both Berends and van de Wal's (2016) and Hinck et al. (2020) algorithms negate the need to prescribe lake boundaries. More sophisticated schemes such as "connected-component schemes" (e.g., Bueler and van Pelt, 2015) are likely to be more efficient for larger-scale computations.
Ice-marginal lakes are commonly hydraulically connected to a glacier and to each other through subglacial drainage systems. Stubblefield et al. (2019) developed and experimentally analyzed a numerical model of the evolution of glacial lakes connected by subglacial channels to analyze filling-drainage cycles of those lakes as the meltwater supply varies. They highlighted that the dynamic and two-way coupling between lake dynamics and glacier behavior affected bed friction, glacier thermal regime and hence ice motion in addition to mechanical buttressing.

Lake Bathymetry and Glacier Bed Topography
The modeling outlined in Meltwater Exchanges strongly depends on topographic data resolution and on lake bathymetry. Indeed, lake bathymetry and glacier bed topography are essential datasets for all numerical models wishing to simulate glacier behavior and lake-glacier interactions. Unfortunately, acquiring lake bathymetry data and bed topography data is very difficult and field data often simply does not exist. The Ice Thickness Models Intercomparison Experiment (ItMix) (Farinotti et al., 2019) has provided a consensus model estimate of ice thickness for all mountain glaciers globally, from which bed topography can be derived. The margins of ice sheets and mountain glaciers are especially poorly constrained because airborne survey instruments struggle to generate reliable data where an ice surface is significantly broken up, the ice substrate is near or at the melting point, and scattering and signal absorption are high due to wet and often debris-rich ice (e.g., Morlighem et al., 2013).

Longitudinal Stress Regime
Ice-marginal lakes fundamentally alter the longitudinal stress and flow regime of a glacier and that effect can propagate kilometres up-ice from the terminus. Field data analysis and remote-sensing of glacier surface velocity by Tsutaki et al. (2019) and by Liu et al. (2020) have both showed that the presence of an ice-marginal lake alters glacier thickness. They both attributed the effects of the lake to be a dynamic thinning mechanism. Tsutaki et al. (2019) additionally conducted some numerical simulations to support their dynamic thinning interpretation but their model had a static lake water level. However, the positive feedback between glacier mass loss, as manifest in ice-margin recession and surface lowering, and meltwater feeding an enlarging lake means that lake levels fluctuate; often by many tens of meters (e.g., Carrivick et al., 2017;Armstrong and Anderson, 2020). Dynamic water levels are technically difficult to achieve and computationally expensive to consider. Similar effects of ice-marginal lakes affecting longitudinal flow extending many kilometres up-ice have been interpreted recently for large lobate outlet glaciers in southern Iceland (e.g., Storrar et al., 2017;Dell et al., 2019;Baurley et al., 2020) and numerically modeled for an ice sheet margin in west Greenland (Price et al., 2008).

Lake Water Temperature
Ice-marginal lake water temperature dramatically affects subaqueous melt rate, but unfortunately is very rarely quantified. There is an unspoken assumption that it remains close to zero due to a predominant influence of subglacial meltwater inputs, but in actuality it can vary considerably between seasons, both in magnitude (to over 10°C) and also in stratification pattern (Haresign and Warren, 2005;Röhl, 2006;Sugiyama et al., 2016;Carrivick et al., 2017;Minowa et al., 2017;Mallalieu et al., 2020;Watson et al., 2020). On a small (meters) scale, thermo-mechanical notch development will be affected by lake water temperature (e.g., Diolaiuti et al., FIGURE 2 | Visualization of outputs from BISICLES ice sheet model (Cornford et al., 2013) used by Sutherland et al. (2020) with starting condition of a glacier in equilibrium (A), then forced by climate only (B) and then by the same climate forcing with the addition of ice-marginal lake effects (C). The modeled lake-terminating glacier has a different terminus geometry and retreats further in the same period of time as the land-terminating glacier.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 577068 2006; Röhl, 2006). On a meso-scale, lake temperature affects subglacial hydrology Stubblefield et al., 2019). On a macro scale (hundreds to thousands of meters) lake water temperature will affect the shape of the floating underside of a glacier and potentially also affect the shape of the entire terminus tongue Sutherland et al., 2020; Figure 2). Subglacial hydrology is affected by lake water temperature at all scales. Lakes can affect local surface energy balance via an albedo effect, as well as a thermal heat capacity effect. Lake water temperature and lake albedo can affect local air movements, and lake water circulation currents can be driven by valley winds. This modification of a local surface energy balance by lakes, and the potential for wind to enhance thermal undercutting and calving, is probably proportional to lake size (Sakai et al., 2009;Sugiyama et al., 2016), but remains unquantified. Very large lakes can act as moisture sources, enhancing precipitation on parts of ice sheets.

Spatio-Temporal Variability and Boundary Migration
Representing the impact of short-lived (decadal) accelerations and decelerations of a glacier in response to ice-marginal lake effects (Sutherland et al., 2020) is difficult within a longer (millennial) glacier evolution model, and not least because those effects might not simply average out. Being able to perform long simulations over large spatial extents must contend with lake size being highly sensitive to the size and position of lake overspill channels, which are on a scale of 10-100 m.
There is a spatio-temporal disconnect between observations made in the field (days to weeks), those from remote-sensing (inter-annual) and numerical modeling concerns (decades to centuries to millennia). This disconnect is a problem because boundary conditions and parameters needed by models vary spatio-temporally. Lake bathymetry, glacier bed topography, lake level and lake water temperature vary in space as well as through time. Multiple ice-marginal lakes will likely have differing elevations, which is quite different to conditions at a marine boundary and makes it difficult to adapt a marine ice sheet model for use with icemarginal lakes (Sutherland et al., 2020). Furthermore, individual lakes will have levels that vary with filling and draining. Convective circulation due to thermal warming of distal water could enhance subaqueous melting. Thermal stratification will likely be an effective contribution to glacier melt seasonally (Sugiyama et al., 2016). The differences between water temperatures at an ice-marginal lake boundary and at a marine boundary are yet to be fully elucidated. A lake-glacier boundary is not a simple interface; it moves in time with lake expansion and glacier retreat, for example. It is complex in that lake water can be hydraulically-connected to a subglacial drainage network.
Calving is assumed to be the dominant effect of a water boundary at an ice-margin. That assumption might be valid for marine-terminating glaciers but remains untested for laketerminating glaciers. Calving models are increasing in number and complexity and they are becoming more and more explicit in their process representation but, because they are orientated toward marine boundaries, they tend to be driven by short-term (intraseasonal) observations. However, some of the water-depth calving laws are derived from data at lake margins (e.g., Benn et al., 2007). The few inter-annual, decadal to millennial scale glacier evolution models are either of a statistical nature; with a stochastic representation of calving, or else they are of a highly-parameterized numerical modeling type. The two most important factors controlling rates and mechanisms of calving at ice-marginal lake boundaries are water depth relative to ice thickness for buoyancy-driven calving (e.g., Boyce et al., 2007;Trüssel et al., 2013) and lake water temperature for melt-undercutting (e.g., Haresign and Warren 2005;Diolaiuti et al., 2006;Röhl, 2006). Buoyancy is different in freshwater compared to marine settings due to water density differences and possibly also ice density differences. A third calving mechanism, which is commonly considered in marine ice sheet models, is hydrofracturing whereby water pressure within surface crevasses causing propagation of that crevasse to the glacier bed, but this calving model has yet to be understood for laketerminating glaciers.
The few studies that have sought to understand glacier or ice sheet evolution with representation of ice-marginal lake effects have adapted marine ice sheet models (e.g., Figure 2). Such an adaptation is not a trivial task, either conceptually for the reasons listed above, but also from mathematical, numerical and computational perspectives. Some of the earlier examples of this sort of work were on the evolution of the Laurentide Ice Sheet and included lakes as a calving term for land-based margins Peltier, 2004, 2006).

Data Available to Validate Models
The previous five sub-sections all concern boundary conditions or data that might be used to parameterize numerical models. Model validation must come from independent quantification of the evolution of lakeglacier boundaries in space and time. For example, a study wishing to numerically model the evolution of the Laurentide Ice Sheet would need geochronological data with good spatial coverage to adequately constrain not only glacier ice margin positions (e.g., Dalton et al., 2020) but also ice-marginal lake shoreline transgressions. Inventories of ice-marginal lake occurrence and areal size are improving in quality and coverage all the time (see Introduction to this article) but looking forwards, there are exciting prospects offered by using cloud computing, such as Google Earth Engine, for efficient automatic identification and geometric delineation of ice-marginal lakes from multiple satellite images and over wide geographical areas.

SUMMARY, CONCLUSIONS AND FUTURE PROSPECTS
This article has emphasized that there is a pressing need to include lake-glacier interactions in numerical models of Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 577068 mountain glacier and ice sheet evolution. It has highlighted that there are a several challenges to overcome for those modeling efforts to be meaningful. The problem of including lake -glacier interactions can be summarized as follows i) determining which processes have significant effects and ii) determining physically-plausible parameter ranges. When moving upscale from single lake-glaciers to whole ice sheets, then also iii) identifying which processes to ignore or simplify, and iv) developing a model framework that can account for and track individual lake size and elevation. This Perspectives paper provides the foundation for future studies aiming to address these problems more thoroughly. Concerted efforts continue by researchers around the World who are applying ice sheet models to ice-marginal lake problems, most especially to reconstructions of the Laurentide Ice Sheet for interest in ice sheet (in)stability and collapse, ice stream development and meltwater pulses. Understanding the last Quaternary (de)glaciation will inform our understanding of the present and future evolution of the Greenland and Antarctic ice sheets. Indeed, we challenge future Greenland and Antarctic ice sheet models to consider lake-glacier interactions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

FUNDING
JLS was supported by a NERC PhD studentship NE/L002574/1. JM was supported by a Graduate Teaching Assistantship from the School of Geography, University of Leeds.