An ice-sheet scale comparison of eskers with modelled subglacial drainage routes

Eskersrecordthesignatureofchannelisedmeltwaterdrainageduringdeglaciationprovidingvitalinformationon the nature and evolution of subglacial drainage. In this paper, we compare the spatial pattern of eskers beneath the former Laurentide Ice Sheet with subglacial drainage routes diagnosed at discrete time intervals from the resultsofanumericalice-sheetmodel.Perhapssurprisingly,weshowthateskerspredominantlyoccurinregions where modelled subglacial water ﬂ ow is low. Eskers and modelled subglacial drainage routes were found to typically match over distances of b 10 km, and most eskers show a better agreement with the routes close to the ice margin just prior to deglaciation. This supports a time-transgressive esker pattern, with formation in short( b 10km)segmentsofconduitclosebehindaretreatingicemargin,andprobablyassociatedwiththin,stag-nant or sluggish ice. Esker-forming conduits were probably dominated by supraglacially fed meltwater inputs. We also show that modelled subglacial drainage routes containing the largest concentrations of meltwater show a close correlation with palaeo-ice stream locations. The paucity of eskers along the terrestrial portion of these palaeo-ice streams and meltwater routes is probably because of the prevalence of distributed drainage and the high erosion potential of fast- ﬂ owing ice. ©


Introduction
Eskers are slightly sinuous ridges composed of glaciofluvial sand and gravel that are deposited in subglacial, englacial, or supraglacial drainage channels (e.g., Banerjee and McDonald, 1975;Brennand, 2000;Storrar et al., 2014a). They can extend for tens to hundreds of kilometres (taking into account small gaps), reach in excess of 50 m in height, and are typically arranged roughly parallel to each other (e.g., Prest et al., 1968;Banerjee and McDonald, 1975;Shilts, 1984;Shreve, 1985a,b;Aylsworth and Shilts, 1989;Clark and Walder, 1994;Boulton et al., 2009;Storrar et al., 2014a,b). Eskers may therefore provide vital information about channelised drainage. This is significant because observations from modern ice sheets reveal the important role of water in lubricating the bed and facilitating rapid ice flow (Bartholomew et al., 2011;Sundal et al., 2011). In particular, the configuration of the subglacial drainage network and how it evolves to accommodate water inputs is critical (e.g., Budd et al., 1979;Alley et al., 1986;Iken and Bindschadler, 1986). Two end-member drainage configurations are typically envisaged (e.g., Walder and Fowler, 1994): (i) an efficient channelised system commonly associated with lower water pressures, lower ice velocities, and higher water discharges; and (ii) an inefficient distributed system (e.g., linked cavities, braided canals, and porous till layer) commonly associated with higher water pressures, higher ice velocities, and lower water discharges. However, because of the difficulties in directly observing the drainage of water at the bed of ice masses, we have a limited understanding of the distribution and geometry of the subglacial drainage network and a lack of data at the spatial and temporal scales necessary to constrain or test subglacial hydrological models (e.g., Hewitt, 2011;Werder et al., 2013).
Presently, investigating eskers located under contemporary ice sheets is physically difficult. Thus the imprint of eskers recorded on the bed of former ice sheets has a clear advantage over data from contemporary ice sheets because we can directly observe the expression of meltwater drainage over large spatial scales. However, despite the use of eskers to reconstruct and constrain ice-retreat histories (e.g., Dyke and Prest, 1987;Margold et al., 2013), very few studies have investigated their pattern at the ice-sheet scale (Aylsworth and Shilts, 1989;Clark and Walder, 1994;Storrar et al., 2014a,b). This is because it is not known whether eskers form synchronously in long conduits (cf. Brennand, 1994) or if they represent a time-integrated signature of drainage deposition throughout deglaciation (e.g., Banerjee and McDonald, 1975;Shilts, 1984;Dyke and Dredge, 1989;Kleman et al., 1997;Hooke and Fastook, 2007). Consequently, the extraction and interpretation of information about where they form in relation to subglacial water generation and ice dynamics is difficult.
A new approach to understanding the pattern of eskers is to compare their distribution and orientation with numerical models of subglacial meltwater drainage at discrete intervals throughout deglaciation. In this study our aim is to compare the expression of eskers on the bed of the former Laurentide Ice Sheet (LIS) with subglacial drainage routes predicted from the results of a numerical ice sheet model. This builds on previous work that has compared eskers and hydraulic gradient routes at glaciers (e.g., Syverson et al., 1994), but represents the first attempt to compare modelled drainage with eskers at the ice-sheet scale and at discrete time intervals.

Mapping esker networks
This paper uses 3749 interpolated eskers (mostly N2 km long) mapped by Storrar et al. (2013). The crestlines of esker ridges were mapped from Landsat 7 Enhanced Thematic Mapper (ETM) + imagery of Canada, which has a resolution of~30 m and~15 m in the panchromatic band. Eskers were typically mapped at a scale of 1:40,000 and were identified based on the criteria set out by Margold and Jansson (2012). Shorter eskers (b 2 km long) were more difficult to identify in Landsat imagery. Comparison with mapping from aerial photographs suggests that~75% of eskers were identified and that 81% of those missed are b2 km long (Storrar et al., 2013).
To enable the effective comparison of eskers with modelled subglacial drainage routes, we used the interpolated esker data set produced by Storrar et al. (2014a). This data set was derived by interpolating a straight line (over short distances in the majority of cases) between aligned esker ridges that appear genetically related (i.e., formed in the same conduit) and merging with the mapped ridges to produce a single esker. It was produced to fill gaps that may have resulted from fragmentary deposition, post-depositional erosion or, submergence beneath lakes and is therefore thought to give a better indication of where the esker-forming conduits were located (Storrar et al., 2014a). We refer to this data set throughout the paper simply as 'eskers'.

Modelling subglacial meltwater drainage
Subglacial meltwater drainage was modelled using the method outlined in Livingstone et al. (2013a,b). Hydraulic potential surfaces (ɸ) of the LIS were calculated from the Shreve equation (Shreve, 1972): where ρ w is the density of water; ρ i is the density of ice; g is the acceleration of gravity; h is the bed elevation; and H is the ice thickness. We calculated the subglacial drainage routes every 500 years for the period between 12 and 7 ka BP, which encompassed the largest retreat distance (hundreds of kilometres) during deglaciation and was over the predominantly hard crystalline bedrock on the Canadian Shield (see Dyke, 2004). The bed elevation data (h) were constructed at 5-km resolution from Gebco_08 digital elevation model (DEM), and the palaeo-ice surfaces and palaeo-bed topographies (corrected for isostasy) were derived from ice-sheet model output from one of the higher probability runs (LT9927) from the ensemble-based analyses of the LIS using the three-dimensional (3D) Glacial Systems Model (GSM) (Tarasov et al., 2012). The GSM includes a 3D thermomechanically coupled shallow ice-sheet model, bed-thermal model, visco-elastic bedrock response, and coupled surface drainage and pro-glacial lake solver. The GSM is calibrated against a large set of observational constraints, including geological and geomorphological evidence and is able to reproduce ice stream locations and ice-margin positions (Stokes and Tarasov, 2010;Tarasov et al., 2012). The LT9927 is from the subensemble of runs used by Livingstone et al. (2013a,b). Their analysis showed that the modelled distribution of subglacial lakes and major drainage routes is a robust result achieved irrespective of the model run used from this subensemble (see Fig. 8 from Livingstone et al., 2013a,b). Given this and the time required for analysing each run, we base our analysis just on LT9927 for this study. The 1°longitude by 0.5°latitude resolution model output was regridded at 5-km cell size. Subglacial drainage routes (i.e., the direction that water flows) were constructed from the hydraulic potential surfaces, using simple GIS routing tools as per Livingstone et al. (2013a,b). Basal meltwater production (cm/y) generated from the GSM was used to weight flow accumulation down subglacial drainage routes. Each cell was given an accumulative basal meltwater value of all the cells that flow into it, hereafter referred to as the 'modelled subglacial flow concentration'. Output cells with a high flow accumulation represent a drainage route along which subglacial meltwater is concentrated.
To allow for basal meltwater production caused by likely subgrid topographic variation, we set a minimum basal meltwater output (0.1 cm/y) in regions of the bed where the temperature was 0 to − 2°C below pressure melting point. Meltwater may also enter the subglacial system from supraglacial sources, although these are not reproduced here because of the difficulty of modelling this process. Thus, we use basal meltwater production simply to indicate where meltwater is likely to concentrate rather than suggesting that all meltwater is necessarily produced at the bed.

Comparison of eskers and modelled subglacial drainage routes
To our knowledge, a comparison between the pattern of eskers and modelled subglacial drainage routes has not been previously undertaken at the ice-sheet scale. Thus, in this analysis, we explore firstorder relationships between the location and orientation of eskers and modelled drainage routes (Fig. 1). Approaches for these two comparisons are described below.

Spatial conformity of eskers and modelled subglacial flow concentration
Output cells with a high subglacial flow concentration indicate regions where large volumes of meltwater are routed, and these represent potential meltwater conduit locations (Fig. 1A). These should correspond to esker locations, as this is where Röthlisberger channels are theorised to form (e.g., Röthlisberger, 1972;Shreve, 1972).
To investigate how the spatial pattern of eskers relates to the routing of subglacial meltwater beneath the LIS, we compared modelled subglacial flow concentration with the esker pattern at 500-year time slices from 12 to 7 ka BP ( Fig. 1A). At each time slice we extracted flow concentration values of all cells that contain eskers and all cells covered by ice in the model domain. We also identified cells that match with the spatial extent of palaeo-ice stream locations recently compiled in Margold et al. (2015) and cells that match the terrestrial portion of the palaeoice sheet bed (i.e., where the mapping was carried out) (Storrar et al., 2013). The probability density function of modelled subglacial flow concentration was calculated for each of the variables extracted and results displayed as a ratio between each probability density function and the probability density function of all the cells covered by ice in the model domain. Statistical significance was evaluated with a binomial test.
To further identify any spatial match between eskers and major subglacial drainage routes we used a flow concentration of N 20 cm/y to identify potential meltwater conduit locations (Fig. 1A). A value of 1 was assigned to cells where the flow concentration exceeded 20 cm/y and a value of 0 to those that did not (see also Livingstone et al., 2013b). This was done for every time slice and the values (0 s and 1 s) then summed together to produce a composite map of potential meltwater conduits and their persistence over time.

Directional conformity of eskers and modelled subglacial drainage routes
Eskers provide a composite history of the former drainage and direction of meltwater flow beneath an ice sheet throughout deglaciation. To try and disentangle what is assumed to be a time-integrated pattern, the match between the orientations of eskers and modelled subglacial drainage routes beneath the LIS were compared at 500-year intervals between 12 and 7 ka BP ( Fig. 1B) (Livingstone et al., 2013a). This comparison was used to investigate (i) how far up-ice from the margin eskers align with the modelled drainage directions; (ii) whether eskers show a better alignment closer to or farther away from the ice margin; and (iii) whether any particular time slice is matched to the entire esker pattern.
A tracking algorithm was developed to compare the directional conformity of eskers and modelled drainage routes (see supplementary information). The tracking proceeds up-ice, moving from the margin of the LIS toward its centre. Tracking was initiated for all eskers at each time slice, with esker-model agreement starting the first time an esker's orientation matched the modelled subglacial drainage direction of an overlying cell and continuing until the agreement first breaks down (Fig. 1B). The length (L) that each esker was tracked and the percentage (f) of overall length tracked ((ΣL tracked /ΣL total ) × 100) was used to investigate the upstream conformity of eskers with modelled drainage directions. We also compared (i) the timing and (ii) the distance from the ice margin between the best and average fits of the modelled subglacial drainage and esker directions at each point along an esker that was tracked in ≥2 time periods. This was used as a measure of when (timing) and where (distance from the ice margin) eskers best fit the modelled drainage directions.
The algorithm has three main parameters that could be varied: (i) the maximum angle of divergence (θ) allowed between the orientations of the esker and modelled drainage direction; (ii) the length scale at which the eskers are smoothed (E); and (iii) the spatial scale (W) used to smooth the modelled subglacial drainage direction grid. Preferred values for the parameters are 45°for θ and a smoothing scale of 5 km for E and W, which is the same resolution at which the model output was regridded and drainage routes constructed. The supplementary section details a sensitivity analysis of θ, E, and W.

Limitations
Our preliminary approach has three major limitations. First, the results use the predictions of a numerical model, which is not reality.  Fig. 2. Spatial relationship of eskers, the domain used to map the eskers (i.e. terrestrial portion of the former LIS bed), and ice streams with varying subglacial flow concentration (ranging from 1 to 100,000 cm/y). Results are displayed as a ratio between each probability density function and the probability density function of all the cells covered by ice in the model domain (i.e., the population). The 20-cm/y drainage cutoff we use to produce Fig. 3 is shown by the vertical line. Values above 1 indicate that subglacial flow concentration happens more often than you would expect by a random distribution. Numbers in brackets are the number of cells used to produce each probability density function. Note how eskers are associated with lower subglacial flow concentrations and ice streams with higher subglacial flow concentrations.
However, the GSM includes basic glaciological physics and is calibrated against a large set of observational constraints (see Tarasov et al., 2012), including relative sea-level data, present-day rates of surface uplift, and an ice-margin chronology (±250-to 1000-year uncertainty) derived from geological and geomorphological evidence (Dyke, 2004). The LT9927 is therefore glaciologically self-consistent and is a sample from the most well-constrained distribution of possible deglaciation chronologies to date. Future extension of this study to a subensemble from the calibration will enable uncertainty quantification of our initial results herein.
Secondly, the difference between the initial resolution of the numerical ice model used to derive the ice surface and bed topographies (~50-km resolution) and the (tens of metre) scale at which eskers were mapped is large. Even the smoothed (5-km resolution) ice surfaces used to generate the hydraulic potential surfaces are relatively coarse, which might artificially create drainage routes that are farther apart than they would be in reality. However, away from the ice margin and ice streams/outlet glaciers, ice sheets are relatively smooth and a 50-km resolution is adequate.
Thirdly, this paper uses a simple representation of subglacial hydrology, which provides a first approximation of where water should seek to flow and pond relative to the ice-sheet model selected. Crucially, the subglacial drainage system and overlying ice are not dynamically coupled, and it does not allow drainage configurations to evolve. We therefore assume that water pressure is equal to ice-overburden pressure uniformly across the bed. This is a reasonable assumption when averaged over coarse ice-sheet scales where geometry is the key control (Livingstone et al., 2013a,b), but unlikely to be valid at smaller scales dominated by local processes (e.g., Gulley et al., 2012). Fig. 2 reveals that the probability of palaeo-ice streams occurring in regions with high subglacial flow concentrations is significantly higher than chance, even if the resolution is reduced ten-fold (binomial test for flow concentrations N20 cm/y, p b 0.001). In contrast, cells that match esker locations show a lower probability for subglacial flow concentration values N100 cm/y. At values b100 cm/y, eskers occur significantly more often than chance (binomial test, p b 0.001 and p b 0.01 if the resolution is reduced ten-fold) and, importantly, display a stronger correlation than the mapping domain. Eskers exhibit a higher probability than ice streams for subglacial flow concentration values b 20 cm/y.

Spatial conformity of eskers and modelled subglacial flow concentration
The spatial pattern between eskers and potential subglacial meltwater conduits (subglacial flow concentration of N 20 cm/y) is shown in Fig. 3. While localised agreement can be identified, for example in the district of Keewatin, the general pattern is of eskers occurring in the gaps between potential meltwater conduits. This arrangement is particularly noticeable in the Northwest Territories and Quebec (Figs. 3B,C). Moreover, the correlation between mapped eskers and the persistence of potential meltwater conduits throughout deglaciation is not clear (Fig. 3). This is demonstrated in the Ungava Bay region (Fig. 3C), where conduits occur repeatedly and in similar locations throughout deglaciation, but eskers are rare.

Directional conformity of eskers and modelled subglacial drainage routes
The length of directional conformity, L, between esker ridges and modelled subglacial drainage routes ranges from b 1 to N 100 km (Fig. 4). The tracking algorithm matches with 65% of the total length of all eskers. The median length of L is 5.4 km, the 5th and 95th percentiles are 0.5 and 42 km, respectively, and the tail of outliers includes matches of N100 km in all except the 7.5-and 7-ka time slices. Length distributions for L are similar at each modelled time slice, but the tail of outliers is shorter compared to the eskers (Fig. 4). Indeed, longer eskers conform less well to the direction of modelled subglacial drainage routes, with the fraction of N~50-km-long eskers tracked b 20% (Fig. 5).
The difference in distance from the ice margin between the best and average fits of the modelled subglacial drainage and esker directions (Fig. 6A) reveals a positively skewed distribution (skew: +3.09), with the peak offset toward negative values (i.e., the fit is better than average when closer to the ice margin). Fig. 6B illustrates the difference in time between the best and average fits of the modelled subglacial drainage and esker directions. The result is a positively skewed distribution (skew: +2.2), with the peak offset toward negative values, i.e., the fit is better than average when eskers are formed during later time slices.
To evaluate whether or not the distribution of L may have been produced by chance, we compared the results against the probability of tracking continuously along a randomly wandering esker (Fig. 5). As θ is ±45°or, 90°of a possible 180°forward arc, the probability p of a 5-km step (i.e., matching E) randomly orientated with respect to the previous one also being tracked is 0.5. The chance of proceeding to the nth step is p n . We ran five stochastic simulations using the length distribution of the eskers. Fig. 5 shows the probability density of the eskers (blue line), the tracked eskers (L) (red line), and an average of the five stochastic simulations (black line). In initial tracking acquisition, and  progressing from the first to the second 5-km step, L does little better than random. This may be because of resolution issues, measurement 'noise', or local factors dominating ice flow at this scale or near the margin. Further tracking (N10 km) occurs much more frequently than chance, which implies an association between eskers and the modelled subglacial drainage direction output.

Spatial association between eskers, palaeo-ice stream locations, and modelled subglacial flow concentration
Our results show that eskers do not coincide with the dominant drainage routes where high concentrations of meltwater (N 100 cm/y) are modelled and routed subglacially (Figs. 2-3). Indeed, there is a highly significant statistical tendency for eskers to occur in regions of low subglacial flow concentration (b100 cm/y) in the gaps between major drainage routes (i.e., left of the vertical line in Fig. 2). The highest modelled subglacial flow concentrations are associated with palaeo-ice stream locations (Fig. 2), where eskers are rare (Fig. 7). In particular, major marine-terminating palaeo-ice stream troughs with large drainage basins, such as Hudson and M'Clure straits, are predicted to have had significant volumes of subglacial meltwater routed down them (Fig. 3). This is unsurprising given that ice streams require warmbased ice, produce large volumes of subglacial water by enhanced frictional melting (Engelhardt and Kamb, 1997;Kamb, 2001), and typically occur in topographic troughs (cf. Winsborrow et al., 2010) where meltwater is focussed. The pattern of meltwater drainage is therefore strongly controlled by the evolution of fast ice flow in the GSM. Significantly, fast ice flow in the model evolves freely as basal ice approaches the pressure melting point, and it is able to reproduce most major palaeo-ice streams (Stokes and Tarasov, 2010) though the shallow ice approximation employed and model grid resolution would caution against overinterpretation of small-scale flow features.
We suggest a number of explanations for the relative absence of eskers on the terrestrial portions of Laurentide palaeo-ice stream beds where subglacial water was concentrated (Fig. 7). First, eskers tend to form on hard beds in Röthlisberger channels, while ice streams are more typically floored by soft sediments that are predicted to drain water through a distributed network (Clark and Walder, 1994). Secondly, increased ice creep associated with high basal velocities will promote the formation of high pressure distributed drainage networks (e.g., Kamb, 1987;Bell, 2008). Thirdly, any eskers that are deposited are unlikely to be preserved because of the enhanced erosional potential of fast-flowing ice (Boulton, 1996).
The formation of eskers in inter-ice-stream regions, characterised by low subglacial flow concentration, may seem counterintuitive because conduit location and growth are thought to be associated with high water discharges (e.g., Shreve, 1972;Schoof, 2010;Hewitt, 2011;Werder et al., 2013). Indeed, it implies that large volumes of subglacial meltwater were not significant in the formation of conduits and eskers beneath the LIS, a point further demonstrated by the absence of suitable palaeo-subglacial lake locations on the relatively flat Canadian Shield (Livingstone et al., 2013a). An alternative source of water is therefore needed, and the most likely option is the penetration of surface meltwater to the bed (e.g., Zwally et al., 2002). Critically, our model does not allow moulin formation or supraglacial inputs, instead assuming that all meltwater was generated subglacially. However, eskers on the Canadian Shield become more frequent later during deglaciation Fig. 6. Histogram showing, for any point that is tracked in ≥2 time periods: (A) the difference in distance from the ice margin between the best and average fits of the modelled subglacial drainage and esker directions. Negative values indicate that the ice margin is closer when the modelled drainage direction best fits the esker direction. Bins are in 50-km intervals. (B) The difference in time between the best and average fits of the modelled subglacial drainage and esker directions. If the best fit between the modelled subglacial drainage and esker directions occurred later than the average fit (i.e., when the ice sheet was smaller) the result is negative. Bins are in 500-year intervals. As θ is ±45°or, 90°of a possible 180°forward arc, the probability p of a 5-km step randomly orientated with respect to the previous one also being tracked is 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (e.g., Prest et al., 1968;Aylsworth and Shilts, 1989), and this coincides with increased rates of atmospheric warming and ice-margin recession (Storrar et al., 2014b). This suggests a genetic link between supraglacial meltwater inputs and esker formation (see also Shilts, 1984;St-Onge, 1984;Aylsworth and Shilts, 1989;Hooke and Fastook, 2007;Burke et al., 2012).
Basal conditions in inter-ice-stream regions may also have been conducive to subglacial conduit and esker formation. For instance, slower ice flow, where creep closure and glacial erosion are reduced, is more suitable for Röthlisberger channel formation and is consistent with observations of eskers in other deglaciated and modern landscapes (e.g., Price, 1966). Moreover, conduits could have been relatively small compared to the size of the esker produced, with growth occurring slowly by the meltout of dirty ice rather than by large discharges through a large channel (e.g., Clark and Walder, 1994;Warren and Ashley, 1994;Hooke, 2005;Hooke and Fastook, 2007). However, sedimentological investigations of eskers comprising horizontally bedded sand and gravel deposits extending all the way across their width suggests conduits can also be the same size as the esker (e.g., Brennand, 1994).

Directional conformity of eskers and modelled subglacial drainage routes
The length of agreement between modelled subglacial drainage directions and eskers is typically b 10 km (Fig. 4) and the ability to track eskers is improved when the ice-margin is closer (Fig. 6A) just prior to deglaciation (Fig. 6B). In addition, no single modelled time slice is able to reproduce the entire esker pattern. This suggests that the majority of eskers were deposited time-transgressively in short segments of conduit close to the retreating ice margin (as frequently inferred; e.g., Banerjee and McDonald, 1975;Shilts, 1984;St-Onge, 1984;Bolduc et al., 1987;Dyke and Dredge, 1989;Hooke and Fastook, 2007). This interpretation is consistent with a supraglacially fed drainage system, where more meltwater is produced and able to penetrate to the bed at lower elevations and in thin ice, near the margin (cf. Hooke and Fastook, 2007;Storrar et al., 2014b).
Agreement between modelled subglacial drainage directions and eskers can also be tracked for hundreds of kilometres up ice in rare cases . These might be explained in two ways: (i) synchronous formation in subglacial conduits that penetrated deep into the ice sheet interior (e.g., Brennand, 2000); or (ii) incremental formation in a stable conduit that migrated up ice during uniform ice retreat. Where the pattern of retreat was uniform, we cannot distinguish between options (i) and (ii). However, we suggest that synchronous formation of very long (N 100 km) eskers is unlikely given the poor correlation with subglacial meltwater concentration (Figs. 2-3) and the high elevation and great thickness of ice in the interior of the LIS (e.g., Peltier, 2004;Tarasov et al., 2012), which would have suppressed the generation and penetration of surface melt to the bed and increased conduit closure by ice creep, respectively.

Further work
This paper uses a simple model to compare the spatial pattern of esker networks with modelled subglacial drainage routes, and as such, a certain degree of fallibility is unavoidable. The next step is to rerun the analysis using a higher model resolution or to take numerical models that can account for dynamic changes in subglacial water pressure and drainage configuration (e.g., Schoof, 2010;Hewitt, 2011;Werder et al., 2013) and apply them to real world situations where the results can be compared to geomorphological data. In particular, Fig. 7. The distribution of inferred palaeo-ice streams (from Margold et al., 2015) and mapped eskers (Storrar et al., 2013) beneath the former Laurentide Ice Sheet. In general, eskers are rarely found on the bed of palaeo-ice streams. sediment erosion, transport, and deposition need to be included in these models to produce the landform imprint that results from subglacial drainage routing and evolution. Our hope is that this preliminary effort offers a useful template for future work.

Conclusions
In this paper, we explore a new approach that compares the expression of eskers beneath the former LIS with subglacial drainage routes predicted from the results of a numerical ice-sheet model at specific time steps during deglaciation.
Eskers are more abundant in non-ice-stream regions where modelled subglacial flow concentration is low. Moreover, we find that the length of agreement between modelled subglacial drainage directions and eskers is typically b 10 km, and the agreement is improved when eskers lie close to the ice margin just prior to deglaciation. This suggests the majority of eskers were deposited time-transgressively in b10-km segments of conduits close to a retreating ice margin and in thin, sluggish or stagnant ice where conduit closure by ice creep is less significant. We suggest that supraglacial meltwater inputs were important for supporting Röthlisberger channel and esker-forming processes.
We find that modelled subglacial drainage routes dominated by high flow concentration are associated with palaeo-ice stream locations. The paucity of eskers along the terrestrial portion of palaeo-ice streams is probably a result of the high velocities and soft-sediment beds that promote distributed rather than channelized drainage, while any eskers that do form are unlikely to be preserved beneath fastflowing ice.