Exploring the ingredients required to successfully model the placement, generation, and evolution of ice streams in the British-Irish Ice Sheet

Ice stream evolution is a major uncertainty in projections of the future of the Greenland and Antarctic Ice sheets. Accurate simulation of ice stream evolution requires an understanding of a number of “ ingredients ” that control the location and behaviour of ice stream ﬂ ow. Here, we test the in ﬂ uence of geothermal heat ﬂ ux, grid resolution, and bed hydrology on simulated ice streaming. The palaeo-record provides snapshots of ice stream evolution, with a particularly well constrained ice sheet being the British-Irish Ice Sheet (BIIS). We implement a new basal sliding scheme coupled with thermo-mechanics into the BISICLES ice sheet model, to simulate the evolution of the BIIS ice streams. We ﬁ nd that the simulated location and spacing of ice streams matches well with the empirical reconstructions of ice stream ﬂ ow in terms of position and direction when simple bed hydrology is included. We show that the new basal sliding scheme allows the accurate simulation for the majority of BIIS ice streams. The extensive empirical record of the BIIS has allowed the testing of model inputs, and has helped demonstrate the skill of the ice sheet model in simulating the evolution of the location, spacing, and migration of ice streams through millennia. Simulated ice streams also prompt new empirical mapping of features indicative of streaming in the North Channel region. Ice sheet model development has allowed accurate simulation of the palaeo record, and allows for improved modelling of future


Introduction
Ice streams dominate discharge in our contemporary ice sheets, routing ice from the ice sheet interior to the margins (Bennett, 2003). Ice streams account for 90% of the discharge of the Antarctic Ice Sheet (Bamber et al., 2000), and a number of ice streams have been identified as key points of vulnerability for the future evolution of the Antarctic (e.g. Favier et al., 2014;Nias et al., 2016;Waibel, 2017) and Greenland (e.g. Joughin et al., 2008;Gillet-Chaulet et al., 2012;Hogg et al., 2016) ice sheets. Given the importance of ice stream dynamics on ice sheet evolution, understanding stream evolution is crucial to elucidating the future of the Antarctic and Greenland Ice Sheets.
Ice stream behaviour within advancing, stable, and retreating ice sheets remains unclear, and accurate simulation of these processes is a key challenge of glaciology (Hindmarsh, 2018 p.625). Observations from contemporary ice sheets and reconstructions of palaeo ice sheets show that ice streams are spatially and temporally variable (Conway et al., 2002;Dowdeswell et al., 2006;Stokes et al., 2009;O Cofaigh et al., 2010;Margold et al., 2015), but the causes of this variability remain poorly understood (Siegert et al., 2004;Horgan and Anandakrishnan, 2006;Peters et al., 2006). Simulations of contemporary ice sheets typically achieve a close fit to observed surface speeds using data assimilation techniques that tune ice sheet parameters such as effective drag (Gillet-Chaulet et al., 2012;Morlighem et al., 2013;Arthern et al., 2015;Gong et al., 2017). While this achieves a close match to observations over decadal timescales, it does not account for long-term changes to bed hydrology, meaning it is expected that such ice sheet simulations will diverge considerably from reality, even given a perfect climate forcing (van der Veen, 1999). And, when modelling regions without present-day ice cover, there are insufficient observations, rendering these data-hungry methods entirely useless. Studies of palaeo ice sheets use a variety of methods to represent basal friction, including relationships with basal elevation (Martin et al., 2011;Seguinot et al., 2016), sediment thickness (Peltier, 2004;Gregoire et al., 2012Gregoire et al., , 2015, or idealised bed classifications (Boulton and Hagdorn, 2006;Gandy et al., 2018). These representations use proxies of bed friction, and the do not capture the evolution of basal hydrology over a glacial cycle.
Improved simulations of ice streams require a representation of the physical mechanisms of ice streaming. Ice streams achieve their enhanced velocity through basal processes, and typically occur over areas of low bed friction. The friction at the bed of an ice stream is likely to evolve due to a number of processes that operate at difference timescales. Over long (millennial) timescales, glacial erosion and deposition can alter the force balance of the bed. Erosion may change the shape of basal obstacles, altering the form drag between the ice and the bed (Schoof, 2002). Sediment deposition may also mask smaller basal obstacles, and change regions where sediment deformation occurs (Bingham et al., 2017;Davies et al., 2018).
The most rapid change to basal friction can occur due to subglacial hydrology, as differences in water routing, subglacial floods, and supra-glacial water inputs can all dramatically alter the volumes and spatial patterns of water at the bed over short (diurnal to decadal) timescale (Nye, 1976;Smith et al., 2007;Hewitt and Fowler, 2008). Detailed kilometre-scale patterning of basal shear stress under ice streams, inferred through observations and inversions, has been deduced for the Antarctic and Greenland ice sheets (Inversion, 2013;Sergienko et al., 2014), who used mathematical modelling to explain them by coupling water and ice flow. Corresponding geomorphological features in the palaeo-record have been observed (Stokes et al., 2016). The thermal regime of the ice, determined by basal friction, geothermal heat fluxes, and englacial ice temperatures, can also influence subglacial hydrology and sliding rates (Payne, 1995). Basal water pressure at the Whillans Ice Stream, West Antarctica, is close to the ice floatation pressure (Engelhardt and Kamb, 1997;Kamb, 2013), and increased surface melting has been linked to increased ice velocity of the Greenland ice sheet (Zwally et al., 2002), both suggesting that subglacial water availability and routing is a key control of ice streaming. Subglacial meltwater can promote ice streaming by saturating basal sediments, and thus encouraging sediment deformation and sliding across bedrock. However, the inaccessibility of basal environments is an obstacle to more robustly describing and understanding subglacial processes (Hewitt, 2011). The representation of ice streaming and subglacial hydrology has been a key area for development for numerical ice sheet models of varying complexity (Greve, 1997;Bougamont et al., 2011;Bueler and van Pelt, 2015). Due to the uncertainty of basal processes simulations usually ignore some processes and idealise others. The palaeo record offers the chance to test whether ignoring and idealising some basal processes is appropriate.
Here, we aim to introduce an idealised representation of subglacial hydrology to simulate the ice stream and ice sheet evolution of the British-Irish Ice Sheet (BIIS). To meet this aim, we address the following objectives: 1) we describe new development of the BISICLES ice sheet model; 2) we employ the BIIS palaeo archive of ice stream locations to test the skill of the ice sheet model in simulating ice stream advance and retreat over the millennial timescale; 3) we discuss the ingredients required to accurately model ice streams of the BIIS.

Ice stream modelling
Accurately simulating ice streams within an ice sheet model is a key concern because of the large influence of ice streams on ice sheet behaviour and stability. Grounded ice streams with marine margins have been identified as points of vulnerability in the future evolution of ice sheets, in particular the West Antarctic Ice Sheet, because ice streams on a retrograde bed slope are prone to Marine Ice Sheet Instability (MISI) (Schoof, 2007). Theoretical investigations of MISI have almost exclusively assumed an ice sheet sliding on bedrock with a viscous power-law relationship between velocity and stress (Hindmarsh and Meur, 2001;Schoof, 2007;Gladstone et al., 2018). However, there is some evidence that a Coulomb friction approach may be applicable close to the grounding line (Iverson et al., 1998;Schoof, 2006). Recent work has shown that a combined Coulomb and viscous power-law approach changes ice sheet profiles, and causes the stable grounding line position to be in shallower water than a plastic power-law only approach (Tsai et al., 2015).
Models of contemporary ice sheets have achieved a close fit to the observed ice surface velocities (Morlighem et al., 2013), using an optimization technique to determine subglacial friction parameters to match surface velocities. In effect, a spatially varying basal friction parameter obtained accounts for unknown bed and englacial properties. Simulations of contemporary Greenland (e.g.  and Antarctica (e.g. Favier et al., 2014;Cornford et al., 2015) are a strong match for the empirical evidence, and may well be a suitable starting point for decadal to centennial future projections, but typically do not allow the basal friction parameters to vary through time. Therefore, centennial to millennial future projections will require ice stream modelling that is capable of ice stream evolution (Aschwanden et al., 2013). Currently, ice sheet models have success modelling ice streams through a theory of thermomechanical instability, where fast ice flow is allowed once the ice has reached the basal pressure melting point (Hindmarsh, 2018). This spontaneous generation of ice streams has been modelled given a flat bed (Payne and Dongelmans, 1997), with regularly spaced ice streams forming due to thermomechanical instabilities of ice sheet flow. Hindmarsh (2009) showed the importance of simulating both horizontal longitudinal and lateral stresses e membrane stresses -when modelling ice stream generation. Idealised experiments that incorporate a representation of subglacial hydrology demonstrate that ice stream formation can be a response to basal water flow rather than ice sheet thermomechanics (Kyrke- Smith et al., 2013).
Applying ice stream modelling to contemporary and palaeo ice sheets is challenging owing to the complexity and scale of real ice sheet beds. The simulation of evolving ice streams has been undertaken in a number of palaeo (e.g Hubbard et al., 2009;Jamieson et al., 2012;Patton et al., 2016;Gandy et al., 2018;Seguinot et al., 2018) and contemporary (Aschwanden et al., 2013(Aschwanden et al., , 2019 studies. Despite progress in simulating ice streams, the skill of models to generate ice streams that evolve over the millennial time-scale has not been adequately tested against the empirical palaeo record.

Palaeo ice streams
Whilst contemporary observations can produce detailed information on current ice sheets, they do not record the multi-decadal, centennial and millennial variability of the ice stream activity and position. Limited contemporary ice stream evolution has been observed, with the most studied example being the change in discharge of ice streams along the Siple Coast of West Antarctica (Retzlaff and Bentley, 1993;Conway et al., 2002). Therefore, the history of contemporary ice sheet observation is not yet long enough to observe the results of significant margin changes, ice temperature evolution, or basal friction evolution. The record from palaeo ice sheets offers an opportunity to test methods for modelling ice streams against an archive spanning millennia, rather than a contemporary snap-shot of ice streaming. The imprint of palaeo ice streams can be determined from a set of well-established diagnostic geomorphological and sedimentological signatures (Stokes and Clark, 1999;Clark and Stokes, 2001;Margold et al., 2015), such as MSGL (Clark, 1993), Trough Mouth Fans (Vorren and Laberg, 1997), and subglacial bedform convergence (Everest et al., 2005). Increased mapping of palaeo-ice stream tracks (Hughes et al., 2014;Margold et al., 2015) offers a record of ice stream and ice sheet evolution over millennia. Efforts to both relatively and absolutely date flowsets of palaeo ice sheets have been extensive (Greenwood and Clark, 2009;Margold et al., 2015Margold et al., , 2018Hughes et al., 2016;Small et al., 2017), resulting in a wealth of palaeo ice sheet mapping.
Whilst the palaeo-record offers a test-bed for ice stream modelling, the modelling of palaeo-ice stream location and evolution can also inform empirical reconstruction efforts. The landform record is increasingly well documented, but will always be an incomplete record of ice sheet behaviour. A model that has sufficient skill to simulate empirically well constrained ice streams may highlight potential locations of unmapped ice streams. In this manner, the model benefits from the extensive testing opportunity provided by the palaeo record, and in turn models can also help target work to extend the empirical record. This creates a symbiotic relationship between numerical modelling and empirical reconstructions.

The British-Irish Ice Sheet
Testing an ice sheet model against the palaeo record requires confidence in the empirical data. The British-Irish Ice Sheet (BIIS) arguably offers the most complete archive of data constraining the behaviour of an ice sheet and several ice streams over millennia. Mapping of glacial features has greatly expanded since the advent of remote sensing techniques (Clark et al., 2004;Smith et al., 2006), complementing significant chronology work (Hughes et al., 2011). Recently, a wealth of offshore data collection has considerably expanded the palaeo archive of the BIIS (e.g Bradwell et al., 2008b;O Cofaigh et al., 2016;Dove et al., 2017). Empirical work spanning decades has cumulatively led to an extensive palaeo archive for the BIIS .
Whilst the record of ice streaming of the BIIS is extensive, it is almost certainly incomplete. The production and preservation of evidence for ice streaming is spatially inconsistent. There is limited evidence of ice streaming on the Irish west coast; despite streaming occurring around the majority of marine margins of contemporary ice sheets. Evidence from the North Sea is also unclear. There, the relatively flat bathymetry should allow ice streams locations not predetermined by topography, allow mobile positioning, and therefore a less clear empirical record. Evidence for similar icestream systems has been found in the Canadian Shield ( O Cofaigh et al., 2010). Holocene erosion and sedimentation during marine transgression in both the North Sea basin and offshore western Ireland is likely another considerable obstacle to mapping landforms of the last glacial cycle.
Numerical simulations of the BIIS have highlighted the importance of ice streaming to the growth and retreat of the ice sheet (Boulton et al., 2003;Hubbard et al., 2009;Patton et al., 2017a). In particular, Hubbard et al. (2009) simulated an ice sheet exhibiting cyclical "binge-purge" behaviour that was controlled by the dynamics of spatially and temporally variable ice streams. Nevertheless, the rich empirical record of the BIIS offers an underused archive to test developments in ice sheet models. Whilst reconstruction-based modelling (Hubbard et al., 2009;Patton et al., 2016;2017a) has proved informative in understanding long-term ice sheet evolution, an idealised approach (Boulton et al., 2003;Gandy et al., 2018), is more suited to an investigation of the processes of ice sheet change.

Methods
We use the BISICLES marine ice sheet model, to which we have added a simple scheme that allows for sliding where basal water is present. BISICLES has previously been successfully applied to simulations of contemporary ice sheets, accurately simulating ice streams using a statistical inversion technique (Favier et al., 2014;Cornford et al., 2015;Nias et al., 2016;Gong et al., 2017). Here, we refer to this new version of the model as BISICLE-S_hydro, which is described below. BISICLES is a vertically integrated ice sheet model with L1L2 physics retained from the full Stokes flow equations (Schoof and Hindmarsh, 2010), including an approximation of membrane stress. These membrane stresses are necessary for producing ice streams of accurate width independent of resolution (Hindmarsh, 2009).
The key addition to BISICLES in this paper is the use of a sliding law which is sensitive to the presence of till water. The sliding law divides the ice sheet into Weertman-and Columb-frcition regions, (Tsai et al., 2015), accommodating both laws by setting the basal shear stress as the minimum of the two stresses, i.e. (1) C is a friction coefficient here a constant equal to 3000 Pa m À1/3 a 1/3 , based on medium values of C from other experiments using BISI-CLES (Favier et al., 2014;Gong et al., 2017;Gandy et al., 2018), adjusted for the flow law exponent. The basal velocity is u b , and m is related to Glen's flow law exponent, n (Glen, 1955;Weertman, 1957). The Coulomb friction coefficient f is~1 (Tsai et al., 2015), and values between 0.33 and 0.5 have been measured for glacial tills in the laboratory (Iverson et al., 1998). We chose as 0.5 as in other studies (Nias et al., 2018), and based on sensitivity with lower and higher values (Fig. S2). The ice pressure is s 0 , and p w is the water pressure. In practice, this means that most of the grounded ice-sheet base experiences Weertman power-law sliding, while a small area near the grounding line will experience Coulomb sliding.
The sliding law for t b is followed when ice is grounded, and t b is 0 once ice is floating. Nias et al. (2018) tested the influence of this sliding law using BISICLES, and found simulations with the Coulomb sliding law experience greater groundling line retreat than simulations with a Weertman sliding law. Ice thermodynamics is considered using an enthalpy transport scheme according to Aschwanden et al. (2012), where an energy density, is conserved rather than temperature T alone; H c is the specific heat capacity, w is water fraction and L is the specific latent heat of fusion. Basal hydrology is approximated according to Van Pelt and Oerlemans (2012), considering balance of water just in a vertical column, ignoring horizontal transport Note, though, that moisture can be transported horizontally within temperate ice. At the base of the ice, frictional heating occurs due to basal sliding. When basal ice is at the pressure-melting point, excess energy is used to melt the ice. The evolution of thickness of the till-stored water layer, W, evolves through a simple equation, where m is the basal melt rate, r w is the density of fresh water, and D is the vertical till-stored water drainage rate, set at 0.005 m/a. The till-stored water drain rate controls the overall water balance and saturation of the till layer, and therefore the magnitude of ice streaming. This value is based on sensitivity experiments (Section S2). Till-stored water does not diffuse horizontally, or have any horizontal drainage, with all of these details approximated by D.
Basal water pressure, p bw , is given by, Here, W 0 is the maximum allowed value of W, set at 2 m, beyond which the till is saturated. A uniform 2 m maximum till-water layer thickness was set, consistent with (Van Pelt and Oerlemans, 2012) and with sensitivity experiments (Fig. S2), allowing simulated ice stream width that best matched empirical data. a is a factor defining the maximum ratio of pore-water pressure (p w ) to overburden pressure, which is achieved in the case of till saturation. We use a ¼ 0.99, in accordance with observations that a~1 (Luthi et al., 2002). g is the acceleration due to gravity (9.81 m s À2 ). When using the Coulomb portion of the sliding law basal shear stress is a function of the basal water pressure, whilst basal water pressure is a function of the amount of basal water present influenced by basal temperature. These adaptations allow for ice stream formation, in a similar way that hydrology has previously been approximated in PISM (Van Pelt and Oerlemans, 2012), with the addition of the combined Weertman-Coulomb sliding law.

Application to the BIIS
We apply the BISICLES_hydro to model an idealised version of the BIIS. The BIIS offers a realistic bed geometry to test BISICLE-S_hydro, along with the extensive empirical record of ice streaming to help test the skill of the model in simulating ice stream positions. These experiments idealise the climate forcing; the experiments are not intended to act as a reconstruction of the BIIS, rather the BIIS acts as a test-bed for BISICLES_hydro to simulate reasonable ice stream width, spacing and position over millennia.

Model-setup
We set up the model domain to cover the entire BIIS (Fig. 2). The easternmost domain edge runs through the North Sea, placed so that as simulated ice reaches the domain edge its normal velocity is set to zero and thus an ice divide forms. This is to represent the effect of confluence between the BIIS and Fennoscandian Ice Sheet in the North Sea during the last glacial cycle (Sejrup et al., 2016;Roberts et al., 2018). This is only an approximation of the effect of confluence, and therefore streaming features close to the eastern margin may be susceptible to domain-edge artefacts. We use an 8 km Â 8 km horizontal grid, with one level of refinement at the ice sheet margin, producing a 4 km Â 4 km horizontal grid. The simulations have 10 vertical levels. The Weertman portion of the sliding law uses a non-linear (m ¼ 3) exponent. We use a crevasse calving model, which models the penetration of basal and surface crevasses (Nick et al., 2010).
To recreate isostatically adjusted bed topography, we adjust modern topography using reconstructions from a Glacio-Isostatic Adjustment (GIA) model (Bradley et al., 2011) (Fig. 2a). GEBCO (Becker et al., 2009) provides modern offshore bathymetry, and SRTM (Farr et al., 2007) provides onshore topography. The Relative Table 1 Summary of experiment set-ups. The results of sensitivity experiments (experiments with the prefix "MAX" or "MIN") are presented in Supplementary Information S2.

Reference
Start point Length (  Sea Level (RSL) change from at 30ka BP is used to deform contemporary topography, maintaining a high-resolution ice sheet bed whilst also accounting for RSL change. RSL is constant through all experiments; the bed does not evolve after the initial RSL correction. Calculation of the surface heat flux uses an idealised elevationconstant surface temperature of 268 K in all experiments. In the majority of experiments, geothermal heat flux is uniformly set at 60 mW m À2 across the domain. In experiment GEOTHERMAL, onshore geothermal heat flux is set to contemporary geothermal heat flux measurements (Busby, 2010;Farrell et al., 2014), which is assumed an appropriate proxy for geothermal heat flux during the last glacial cycle. Owing to the absence of sufficient measurements, offshore geothermal heat flux is set at 60 mW m À2 .
The simulations aim to grow and shrink an ice sheet in an idealised cycle, rather than provide a reconstruction of the BIIS. Although a more idealised surface mass balance forcing can be used, the scheme must force a simulated ice sheet that resembles the actual evolution of the BIIS well enough so as not to impede comparisons of modelled and empirically reconstructed ice streams, as ice stream position is influenced by ice sheet geometry. The simulations are split into three subsequent experiments with varying surface mass balance forcing (Table 1). First, the ice sheet is grown for 20,000 years from no ice using an ice mask (Fig. 2b). A surface mass balance of 0.3 m/y is applied in the ice mask, and a mass balance of À5.0 m/y applied outside the ice mask. A 2-phase spin-up is used, with the hydrology adaptations discussed in section 2, introduced at model year 10,000 to preserve computer time during the original build-up of the ice dome. From the end of the SPIN-UP experiment, an idealised climate is used to advance and retreat the ice sheet. During the ADVANCE phase, the surface mass balance is defined as; where m is the annual surface mass balance (m/y), and h is the surface elevation of ice. This forcing allows an advance that reaches a maximum extent after 10,000 years which is comparable to empirical reconstructions (Clark et al., 2012). Every 25 years through the ADVANCE experiment, a mask is created based on the modelled output (Fig. 3). The mask sets values of 0.3 where the ice sheet elevation is > 500 m, 0.0 where the ice sheet elevation is 500 m, and À2.0 where the no ice is present. Then the retreat is forced using the masked SMB maps in reverse order, i.e. the mask produced using the 9,000 model year extent forces the 11,000 SMB, the 8,000 year mask forces the 12,000 SMB etc. (Fig. 3). The 500 m ELA provides a balance between too slow and too rapid retreat of the ice sheet in the retreat phase. The surface mass balance forcing is updated every 25 model years. This forcing method creates a near-symmetric pattern of advance and retreat.

Model-data comparison
We compare modelled and empirically reconstructed ice streams, considering ice stream position, spacing, flow direction, and evolution. To allow for consistent identification of ice streams in the simulations we define an ice stream as a region with surface velocity exceeding 500 m/y 8 km from the margin. Qualitative comparisons between modelled and empirical ice streams were made based on position and spacing, and flow directions. This was completed for topographically well-focussed ice streams since these are the locations where the flow direction can be empirically reconstructed with most certainty. The empirically reconstructed palaeo flow direction is determined as the mean orientation of the topographic trough centreline. Together, these were used to create a map of empirically reconstructed flow directions for known ice streams of the BIIS.
We use the tool developed by Li et al. (2007), Automated Flow Direction Analysis (AFDA), to compare empirical and model derived flow directions. AFDA calculates the mean residual angle and variance of offset between modelled and empirically derived iceflow directions. A threshold of <10 mean resultant vector and <0.03 mean resultant variance was used to assess model-data agreement, derived from values previously used in the literature Ely et al., 2019).

Model results
The skill of the simulations to produce ice streams of a reasonable width, position, and flow direction is highly dependent on the "ingredients" included in the simulations. A simulation with no basal hydrology coupling does not simulate ice streaming (Fig. 4e). Simulations using BISICLES_hydro achieve ice stream generation in manner that is a good match to empirical data, and demonstrates good consistency in the simulated ice stream position, spacing and width across a range of horizontal resolutions (Fig. 4aed). Consistent ice stream width means that ice catchments are also consistent across varying horizontal resolutions, allowing margin evolution to not vary considerably between resolutions (Fig. 4aed.). The minimal resolution dependency could be explained by either the inclusion of membrane stresses (Hindmarsh, 2009), or topographical constraints on both ice stream position and width. With increasing resolution there is convergence of simulated velocity (Fig. 4f). At 4 km resolution and beyond the pattern of ice streaming is qualitatively and quantitatively (Table S4) similar.

ADVANCE and RETREAT
At the end of the SPIN-UP experiment, there are four major ice streams draining the ice sheet (Fig. 5a). As the ice sheet advances during the ADVANCE experiment, the position of some ice streams migrate, and the number of ice streams increases (Fig. 5aee). Some ice streams remain spatially and temporally consistent, concurrent with the position of ice streams marked in Fig. 1 At the maximum extent of the ice sheet, at 10,000 model years (Fig. 5e), the number of ice streams around the ice sheet margin has increased during the ADVANCE phase. Ten ice streams have formed on the western-board of the ice sheet, which has a marine margin at the Atlantic Ocean. The areal extent here is controlled by the extent of the continental shelf, as ice cannot extend beyond the shelf break. A similar effect has been modelled on the southern margin of the Laurentide Ice Sheet, where advance into deeper lakes initiates increased calving (Cutler et al., 2001). Ice streams are less prominent along the southern margin, and migrate laterally during the ADVANCE phase.
During the RETREAT experiment (model years 10,000e20,000) all the non-transitory ice streams of advance remain consistent spatially and temporally. Streaming along the southern margin increases, and a large ice stream forms in the Irish Sea (Fig. 5f). Smaller ice streams form onshore, concurrent with the Vale of Cheshire and York, although these stop streaming by 15,000 model years ( Fig. 5g and h). At 12,500 model years, two large ice streams drain the ice sheet on the eastern domain edge in the North Sea. These large features are an artefact of the ice sheet being separated from the domain edge following confluence during the ADVANCE experiments, and therefore are not comparable with the empirical record. Once the ice sheet is separated from the eastern domain edge (~15,000 model years), ice streams in the North Sea area have a comparable size and spacing to other ice streams simulated around the margin.
During the ice sheet advance, steeper margin and shallow interior surface slopes occur along the southern margin (Fig. 6), as is expected on an advancing margin. In the retreat phase the SMB imposes melting at the southern margin of the ice sheet. This lowering steepens the ice surface at the ice stream onset zone, accelerating ice streaming and resulting in a shallower surface profile downstream of the onset zone. This means that during the advance phase there is only extensive streaming along the northern margin, whilst streaming also occurs along the southern margin during the retreat phase (Fig. 6). The feedback between SMB and ice stream behaviour has been demonstrated numerically (Robel and Tziperman, 2016), and this effect has been shown to promote rapid deglaciation. Experiments of the evolution of the BIIS with a realistic climate forcing also produce ice streams with high temporal variability (Hubbard et al., 2009).

Variable geothermal heat flux
Introducing a spatially variable geothermal heat flux produces limited changes to the pattern of ice streaming simulated for most of the ice sheet. The width and spacing of the ice streams is almost identical (Fig. 7). The dominant ice streams identified from the ADVANCE and RETREAT experiments (i.e. coincident with the Minch Ice Stream, and the Barra Fan Ice Streams) remain constant in the variable geothermal heat flux experiments. The region of highest geothermal heat flux in the domain, in South-west England, is not Fig. 5. Ice velocity at 2,500 model year intervals during the ADVANCE (aee) and RETREAT (feh) experiments. Orange boxes highlight ice streams with significant empirical evidence (Fig. 9). Green boxes highlight ice streams with empirical evidence presented in this research (Figs. 11 and 12). The green contour shows the 0 m bed elevation, and the ice sheet grounding line. An animation of the experiment is contained in the supplementary information (S3). covered by ice, both in these simulations and according to empirical reconstructions (Clark et al., 2012), so has no impact on the ice stream dynamics simulated. Other comparative "hot spots" also produce a minimal effect, despite being ice covered. The higher geothermal heat flux over the North of Ireland does not cause significant change in the margin position as the ice sheet advances over Ireland (Fig. 7b,f).
The primary change in the ice stream configuration is on the west coast of Ireland, where the geothermal heat flux is comparatively high (Fig. 2c). In the spatially constant geothermal heat flux experiment an ice stream is simulated coincident with Donegal Bay (Fig. 1[7]), with a second ice stream to the north in the Sheephaven ice stream region (location indicated in Fig. 1[6]). In the spatially variable geothermal heat flux experiment, an ice stream is still simulated in Donegal Bay, but to the south in Clew Bay rather than the north (Fig. 1[8]). In this region of the ice sheet, there are no significant bathymetric troughs, suggesting that variable geothermal heat flux becomes a comparatively more significant control on ice stream form and location in western Ireland than the rest of the BIIS. This demonstrates that geothermal heat flux can be an important control on ice stream location, but only in regions where other more significant controls, like topographic troughs, are not present. Including a spatially variable geothermal heat flux is not a necessary ingredient for modelling the majority of ice streams of the BIIS when using the physics in BISICLES_hydro, but is more influential on the western Irish coast.

Model-data comparison
We compare simulated and empirically-based ice streams using both qualitative comparisons of position and spacing and a quantitative tool. Ice streams of the BIIS are split into three categories; ice streams that are simulated by the model and have empirical evidence (group A), ice streams that are simulated by the model but have no empirical evidence yet published (group B), ice streams that are not simulated but have been reconstructed empirically (group C). The majority of ice streams simulated here (11/17) are supported by empirical evidence (Table 2).

Group A
Group A are the eleven BIIS ice streams simulated here and supported by the empirical record ( Table 2). Four of the most consistent ice streams in the simulations coincide with the Fig. 6. Transect of the ice surface at a period before and after the onset of significant ice streaming in the Celtic Sea (10,925 and 11,100 model years respectively). Grey curves show the ice surface at 25 year intervals between the two snapshots.
locations of the four ice streams with the strongest empirical evidence. Each of the ice streams (the Barra Fan Ice Streams, the Minch Ice Stream, the Rona Ice Stream, and the Foula Ice Stream) are associated with a clear bathymetric trough (Fig. 9), and terminate at major Trough Mouth Fans (Bradwell et al., 2008b). These are thick accumulations of sediment fed across an ice stream marine margin (Vorren and Laberg, 1997). The bathymetric troughs and associated Trough Mouth Fans are strong empirical evidence of sustained ice streaming, which the model here shows significant skill in simulating.
A number of ice streams in Group A form in bathymetric troughs, and therefore locational (or colocation) between the reconstructed and simulated ice streams was likely. However, some of the Group A ice streams form in relatively minor bathymetric troughs, like the North Sea Lobe, the Forth Ice Stream, and the Moray Ice Stream. In these examples, the simulations show significant skill in simulating empirically reconstructed ice streams with subtle bathymetric confinement. Although bathymetry does not directly control the position of these ice streams, neighbouring trough controlled ice streams can influence their position. This means that the skill of the model in simulating ice streams in significant bathymetric troughs allows the model to show skill in simulating ice streams away from bathymetric troughs.
The empirically well-constrained ice streams in bathymetric troughs allow a quantitative comparison between simulated ice streams and empirically reconstructed ice stream paths using AFDA (Li et al., 2007;Napieralski et al., 2007;Ely et al., 2019). The Automated Flow Direction Analysis (AFDA) tool calculates the mean residual vector and variance between simulated and reconstructed ice flow. Fig. 10 shows the resulting mean residual vector and variance of the flow direction for the duration of each ice stream's simulation. The Rona Ice Stream (Fig. 10b) is the strongest match between simulated and empirical evidence, with an average mean residual vector of 5 , and only brief excursions above the datamodel match criteria. The Minch and the Barra Fan Ice Streams also perform reasonably well during periods of the ADVANCE and RETREAT experiments.
The Foula ice stream is the shortest-lived ice stream, streaming from 8,500 to 15,000 model years ( Fig. 5e and f), and matches empirical flow directions at the start and end of the ice stream's existence, but there is no model-data match during the majority of the ice stream's occupancy. There is significant flow direction change of the Foula ice stream during the simulation, as the shape of the margin changes during confluence with the eastern domain edge.
The results from AFDA confirm the qualitative conclusion that for many ice streams there is often a strong match between model and empirical data. The empirically best-constrained ice streams of the ice sheet form in significant bathymetric troughs (Fig. 9), which the model demonstrates significant skill in simulating. Whilst there may be a weaker model match for ice streams that are less bathymetrically controlled, like the North Sea and western Irish margin, the empirical reconstruction of flow direction is also less well constrained in these regions. The AFDA scores are sensitive to variations in the empirical reconstruction, so poorly constrained ice streams are not suitable for meaningful comparison.
A number of simulated ice streams are more prominent during the retreat phase than the advance phase. These ice streams include the Moray Firth Ice Stream, the North Sea Lobe, and the Irish Sea Ice Stream. The retreat forcing increases the propensity of the ice sheet to stream along the southern margin (Robel and Tziperman, 2016). The simulated ice streams of the retreat phase generally form in positions with strong empirical evidence, like the North Sea Lobe (Dove et al., 2017). Ice streams that are simulated only during the retreat of the ice sheet form in less significant bathymetric troughs that the Minch, Barra Fan, Foula, and Rona Ice Streams. It seems, therefore, that these have a lower propensity to stream, until triggered by the ice sheet retreat.

Group B
Simulated ice streams without published empirical evidence form the second largest group of ice streams of the BIIS; group B (Table 2). Empirical evidence for ice streaming is inconsistently preserved by subsequent erosion and deposition. Because of this, the empirical record is, and will always be, incomplete. Given the apparent skill of the model (Group A), that ice streams may be simulated which are currently unmapped, is not a basis for discounting the model results.
The presence of simulated ice streams with no published empirical evidence stimulated us to conduct a re-examination of the empirical record in these areas. This is most prominent for the Irish west coast, where a number of ice streams are simulated, but without empirical evidence for streaming reported. Here, high rates of sedimentation during the Holocene may be obscuring bathymetric troughs and subglacial lineations. In the future, reflection seismic data may reveal these features.
In other areas, improved bathymetric data coverage allows new mapping to reveal evidence of ice streaming in regions simulated here. A North Channel Ice Stream, flowing in the Irish Sea west of the Isle of Man, is simulated in the advance and retreat of the ice sheet (Fig. 5). Recent bathymetric multibeam data, from the UK Hydrographic Office with point data gridded into an 8 m horizontal grid, allows major subglacial lineations in this region to be mapped (Fig. 11). Two hill-shades were created for the gridded bathymetric data, illuminating from 045 to 315 to avoid azimuth biasing (Smith and Clark, 2005). Hill-shaded bathymetry is overlain with translucent bathymetry to aid with geomorphological feature mapping. The mapped subglacial lineations are elongate (typically 2e6 km long and 100s of m wide), indicative of high ice velocities (Clark, 1993;Stokes and Clark, 1999;Spagnolo et al., 2014;Ely et al., 2016). This example demonstrates how model results can motivate and steer empirical work, here leading to the identification of a North Channel palaeo-ice stream that has not been previously documented.
A St Kilda ice stream (flowing just south of the Scottish archipelago St Kilda) is also consistently simulated throughout the glacial cycle. Here, the bathymetric data does not provide evidence for subglacial lineations, but there is some circumstantial evidence for ice streaming (Fig. 12). A small bathymetric trough is evident, encouraging ice flow between South Uist and Barra. The region is predominantly exposed bedrock, potentially caused by the Fig. 8. Ice velocity at 10,000, 12,500, and 15,000 model years. Locations of key ice streams empirically reconstructed in the literature are highlighted with black arrows, with ice streams numbered as in Fig. 1 and Table 2. The wide ice streams in the North Sea are caused by a domain-edge effect. extensive ice stream erosion (e.g. Bradwell et al., 2008a;Newton et al., 2018). Whilst the evidence for the North Channel Ice Stream is direct, the empirical evidence for a St Kilda Ice Stream remains circumstantial. However, the model results could encourage more targeted empirical data collection to investigate the evidence for a palaeo-ice stream.
Overall, whilst the majority of ice streams simulated have strong evidence in the empirical record, there are a handful of ice streams simulated with no empirical evidence. This offers an opportunity for targeted empirical data collection and analysis to find evidence for or against ice streams simulated in these experiments, which is supported by increasing quality of bathymetric products (Becker et al., 2009;Calewaert et al., 2016), and greater collection of shallow reflection seismic data (O'Brien et al., 2016;Stewart, 2016).

Group C
A small number of BIIS ice streams have empirical evidence but   are not simulated here (Group C). This includes the Tweed Ice Stream (Fig. 1[18]), evidenced by extensive mapping of streamlined subglacial bedforms (Everest et al., 2005). The flow of the Tweed Ice Stream is dependent on margin positions in a complex region of the ice sheet (Livingstone et al., 2015). The simulations here are not intended to be a margin reconstruction of the BIIS, and do not include the complex margin history expected in this region.
The Celtic Sea Ice Stream is the ice stream of the BIIS with the greatest width that is not consistently present in the simulations. The southern margin of the simulated ice sheet in the Celtic Sea never reaches as far south as empirical reconstructions (Scourse et al., 2009). The limited model extent compared to empirical data in this region may be due to a lack of persistent ice streaming during the advance phase, the lack of a simulated large surge event, or because the idealised SMB forcing in this case is sufficiently different than reality in this sector. Whilst the empirical evidence is able to constrain the extent of the Celtic Sea Ice Stream (Scourse et al., 2009;Praeg et al., 2015), the flow geometry is difficult to reconstruct owing to a lack of subglacial bedforms. The differences in the empirical evidence between the Celtic Sea Ice Stream and other ice streams of the BIIS, along with the low model skill compared to high skill for other ice streams, suggests that the Celtic Sea Ice Stream was mechanistically different to the other ice streams of the BIIS.

Required ingredients for modelling ice streams
Here we discuss which model ingredients are crucial to successfully model ice streams, which are desirable, and which can be deemed as not important. The most important model ingredient in these experiments is the representation of idealised subglacial hydrology, in these experiments as a till-water layer coupled with the Coulomb portion of the sliding law, as without this the ice sheet model does not spontaneously generate ice streams (Fig. 4e). An adequate horizontal resolution is also an important model ingredient, the Minch Ice Stream modelled down to 4-1 km resolution (Fig. 4d) achieves an improved AFDA score than when it is simulated at 8 km resolution (Table S4). However, lower-order models experience more significant resolution dependency (Hindmarsh, 2009), and if an experiment was also considering grounding line dynamics using a model of sufficient physical complexity, like BISICLES, could also be considered a crucial model ingredient.
Depending on the context of the experiments, additional model ingredients are required. The representation of SMB has an influence on ice stream behaviour, evident by the increased streaming along the BIIS southern margin when SMB became negative (Fig. 6), and previous experiments (Robel and Tziperman, 2016). An accurate SMB also has a number of secondary effects. For example, ice streams of central northern England were likely not simulated here because they are dependent on complex margin and ice dome changes that cannot be simulated with an idealised SMB scheme. A spatially variable geothermal heat flux was also determined to be non-crucial for the BIIS but might be important for simulating ice streams of other ice sheets, like the Northwest Greenland Ice Stream (Rysgaard et al., 2018). Previous simulations of ice streams of the BIIS also used a spatially uniform geothermal heat flux (Hubbard et al., 2009). Finally, some model ingredients were either ignored or highly idealised in this study and are not important in this experimental context. For example, a good match to empirical data was achieved despite a spatially uniform Weertman coefficient, Coulomb friction angle, maximum till-water depth, and till-water drain rate. Varying these parameters does change the pattern of ice streaming (Fig. S2), but representing these factors as spatially variable is not a necessary model ingredient in this case. However, a more realistic representation of the subglacial environment, including factors such as Fig. 13. a) Ice Velocity at t ¼ 4,000 years. b) Till water thickness at t ¼ 4,000 years. variable bed geology, and/or improved process understanding may prove to be a more influential model ingredient for other ice sheets, and could help improve the simulation of the Celtic Sea Ice Stream.
The horizontal transport of meltwater is a process that remains unrepresented in the model. The supraglacial, englacial, and subglacial transport has been observed to redistribute water on short timescales (Andrews et al., 2014), and has been identified as a key control on ice sheet velocity (Zwally et al., 2002;Hewitt, 2013). Meltwater transport also evolves annually (Chandler et al., 2013), and is expected to be a mechanism to cause ice stream evolution. Although the experiments here achieve a close fit to empirical data without a representation of horizontal meltwater transport, future model development should include this mechanism and test for its influence in a variety of ice sheet contexts. Winsborrow et al. (2010) reviewed the literature on palaeo and contemporary ice streams, and compiled seven factors that control ice stream location. They are, in the order of importance that Winsborrow et al. (2010) proposes, topographic troughs, marine margins, soft beds, abundant meltwater, smooth beds, high geothermal heat flux, and topographic steps. Whilst ice streams will spontaneously form due to thermo-mechanical coupling of ice flow (Payne and Dongelmans, 1997), these seven factors controls the relative position of ice streams. Of the seven factors, only soft beds are not represented in these simulations as bed friction coefficient is idealised to be uniform across the domain.

Ice stream controls
For the BIIS, the influence of topographic troughs and steps are considered together, as the features co-exist at the ice sheet bed. Supporting the Winsborrow et al. (2010) hierarchy of controls, these simulations suggest that bed topography is the primary control of ice stream location for the BIIS. The primary large, and spatially and temporally consistent ice streams in the simulations, are all located in bathymetric troughs; the Minch Ice Stream, the Barra Fan Ice Streams, the Rona Ice Stream, and the Foula Ice Stream (Fig. 9). However, not all simulated ice streams form in well-defined topographic troughs, like the North Sea Lobe, Forth Ice Stream, and the ice streams of the Irish west coast. Ice streams without welldefined topographic control are more likely to be sensitive to weaker controls on ice stream location, demonstrated by the change in ice stream location in Northwest Ireland when considering a varying geothermal heat flux (Fig. 7). Ice streams outside a topographic trough also exhibit greater resolution dependency (Fig. 4g). An idealised experiment of a square ice sheet on a flat bed (Fig. S1) also shows greater resolution dependency than on a real bed (Fig. 4), suggesting that these ice streams may also be more sensitive to resolution.
The position of topographically controlled ice streams has an indirect control on ice streams with weaker topographic control. It has been shown here (Fig. S1), and in previous research (Payne and Dongelmans, 1997;Hindmarsh, 2009) that even on a flat bed regularly spaced ice streams will form around the margin of an ice sheet. Therefore, whilst topographic troughs will determine the location of some of the ice streams of the BIIS, other ice streams would be expected to form between these ice streams even on flat beds. Ice sheet models on flat beds predict regular spacing of ice streams; but on a real bed with topographic variation, some ice streams would be anchored in topographic troughs and the predominant spacing control then fixes the position of other nontopographic ice streams.
The simulations also support the expected strong influence of marine margins, with the most spatially and temporally consistent ice streams occurring along the western board of the ice sheet, with a marine margin in the Atlantic Ocean. Many simulated ice streams with a marine margin are also in bathymetric troughs, although the ice streams of the western Ireland coast have a marine margin without large troughs. Calving at the margin of the ice sheet lowers the surface profile, allowing the capture of more ice from the catchment, and thus increased ice streaming. Ice streams that advance and retreat onshore, for example through central and southern England (Fig. 5d and e), do not stream constantly, unlike the majority of ice streams with a marine margin. Therefore, there is evidence from these simulations that a marine margin is a strong control of behaviour of BIIS ice streams, but not necessarily ice stream location.
The saturation of the till-water layer is a control on the water pressure, which is in turn a variable in the Coulomb portion of the sliding law. The saturation of the till-water layer is influenced by the geothermal heat flux, maximum till-water layer thickness, and the till-water drain factor. Whilst there is a strong relation between the location of modelled ice streams and regions of a saturated tillwater layer (Fig. 13), the till-water layer saturation is not the sole control on ice stream position. On the southern and eastern margin there are a number of regions of saturation, whilst ice streaming does not occur. However, no ice streams are apparent in regions with low till-water layer saturation. The increased velocity of an ice stream increases frictional heating, and therefore increases the saturation of the till water layer. The relationship between the saturation of the till-water layer and ice velocity means evidence of cold-based ice (Bierman et al., 2015;MacGregor et al., 2016) can be a strong constraint on model results.
The relatively limited influence of the geothermal heat flux is highlighted by the experiments ADVANCE_GEOTHERMAL and RETREAT_GEOTHERMAL, suggesting that the spatial distribution of geothermal heat flux is only a weak control on the formation and location of ice streams. The most significant difference to ice streaming made by a spatially variable geothermal heat flux is in western Ireland, where the influence of bathymetric troughs is weak. The relatively small influence of geothermal heat flux has been supported by Winsborrow et al. (2010). However, it is important to note that there is only a limited range of geothermal heat flux values across the domain. Geothermal heat flux could be an important control of ice stream formation and location in domains with weaker other controls on ice stream location and a wider range of geothermal heat flux values. For example, it has been suggested that the location and unusual geometry of the Northeast Greenland Ice Stream may be influenced by a geothermal heat flux hot spot at the ice sheet bed (Rysgaard et al., 2018).

Model evaluation
The comparison between the simulated and empirically reconstructed ice streams can be used to evaluate the skill of the ice sheet model at producing ice streams in reasonable locations. The qualitative and quantitative similarities between the simulated ice streams and the ice streams recorded in the empirical data suggests that the model has significant skill. The simulations consistently produce ice streams in areas of strong empirical evidence. Therefore it would be reasonable to apply this to other palaeo ice sheets where the empirical record is less complete. If the model can consistently simulate expected ice stream and ice sheet evolution over millennia for a diverse range of palaeo ice sheets, the model would be expected to have sufficient skill to project medium-and long-term evolution of the Greenland and Antarctic Ice Sheets.
However, some mechanisms not represented in the model, which could be influential to ice stream formation and location, including water drainage from the ice surface to the bed (Das et al., 2008;Krawczynski et al., 2009), and ice stream freeze on and shutdown (Christoffersen and Tulaczyk, 2003). These simulations include no representation of spatially variable maximum till-water layer thickness, deformability, or porosity. We consider these geological controls on ice stream formation and location as second order factors given the similarities between simulated ice streams and empirical ice streams show here, but this will not necessarily be the case for all ice streams in all ice sheets. For example, the poor representation of the Celtic Sea Ice Stream in these simulations may be the result of the lack of representation of key geological parameters that may help produce large surging lobes. For example, the lack of horizontal meltwater transport in the simulations still allows good model-data fit for the majority of ice streams, but may be a more important process for the Celtic Sea Ice Stream. A thermomechanical induced surge would also be highly sensitive to model set-up and subsequent evolution, and would be a challenge to reproduce in simulations. Overall, the method described of modelling palaeo-ice streams demonstrates significant skill in producing expected ice streams, especially given the lack of spatially variability in geological parameters in our experiments.

Model-data symbiosis
It is important that palaeo-ice sheet reconstructions use the wealth of information available from both the empirical record and numerical modelling experiments (Andrews, 1982), but in reality this has proved difficult to implement. Comparison between simulated and empirical data is generally limited, with the exception of the comparison between ice sheets and RSL using GIA modelling (Kuchar et al., 2012;Patton et al., 2017b). Despite this, tools to aid the comparison of model and empirical data have been developed (Li et al., 2007;Napieralski et al., 2007;Ely et al., 2019). As simulated and empirical data becomes more extensive, the foundations are in place for a symbiotic relationship between simulated and empirical data, where models benefit, and benefit from, empirical data, resulting in improved palaeo reconstructions, and future ice sheet projections.
Here, the empirical palaeo archive has proved to be a good test of the skill of an ice sheet model's ability to produce ice streams. Empirical data permit qualitative and quantitative assessment of ice sheet model skill at simulating ice streams of reasonable position and spacing. Simulations with a more advanced forcing, primarily a less idealised climate forcing, can benefit from a more extensive suite of model-data comparison tools than the single tool used here. For example, Ely et al. (2019) describe a suite of three tools to compare simulated palaeo ice streams with empirical evidence, which compare margin position, flow direction, and advance and retreat timing.
Simulated data can also complement empirical data collection and analysis by indicating the potential locations of palaeo ice streams. An intriguing candidate, highlighted here, is the St Kilda ice stream (Fig. 12), which is simulated consistently throughout the advance and retreat of the ice sheet but empirical evidence for this ice stream remains circumstantial. The simulations also point to the potential for west coast of Ireland ice streams, with potential empirical evidence subsequently covered in significant Holocene deposits. This means that surficial evidence for these potential ice streams may not exist, but could be mapped through extensive reflection seismic data (Graham et al., 2007;Knutz et al., 2019;Emery et al., 2019), which would be important for the dynamics of a significant marine portion of the BIIS.
Modelling results can also assist reconstructions of ice sheet behaviour in regions where there is limited empirical data. For example, despite extensive empirical efforts, the deglaciation dynamics of the North Sea remain uncertain (Roberts et al., 2007;Sejrup et al., 2016;Dove et al., 2017). An ice sheet model with sufficient skill to model evolving ice stream dynamics could help test the feasibility of empirical reconstructions. This process of an ice sheet model being used to identify physically plausible glaciation patterns could also be applied to other empirically uncertain regions of the ice sheet, like the Porcupine Bank , or the Celtic Sea Ice Stream (Scourse et al., 2009).

Contemporary ice stream modelling
The evolution of ice streams is a considerable uncertainty in projections of Greenland and Antarctic response to climate change (Vaughan and Arthern, 2007), and thus is a considerable uncertainty in projecting long-term global sea level rise (Bamber and Aspinall, 2013). The method of inverting ice stream parameters from surface velocities, which is common in studies projecting future ice sheet evolution (Morlighem et al., 2013), is not suitable for medium-to long-term projections because it does not properly account for the evolution of bed hydrology as the ice sheet evolves (Aschwanden et al., 2013). A statistical inversion of the surface velocity of an ice sheet will produce a close match to the observed ice sheet in the short-term. However, over centennial timescales the observed and modelled ice sheets would diverge because of the evolving nature of ice streams, even given a perfect climate forcing.
Therefore, to project the medium-to long-term evolution of the Greenland and Antarctic Ice Sheets, an ice sheet model is required that simulates ice stream location, spacing, and evolution accurately without inverting ice stream parameters from contemporary ice sheet surface velocity. BISICLES_hydro simulated ice streams without requiring inversion, and accurately simulates ice stream position, spacing, and evolution as suggested by empirical data from the palaeo record. Future work could validate ice velocities simulated with this scheme against measured contemporary ice velocities of the Greenland and Antarctic Ice Sheet. It would also be useful to test the model on sectors of palaeo ice sheets with a less dominant topographic control on ice stream locations, like the southern margin of the Laurentide and Fennoscanidan ice sheets, and to experiment with a more realistic transient climate and ocean forcing. These sectors would also provide a stronger test of the model's ability to simulate terrestrial ice streams, which are infrequent for the BIIS. Assuming sufficient skill, this model would have the potential to be used to project the medium-to long-term evolution of the Greenland and Antarctic Ice Sheets.
Further development could also be made in the representation of englacial and subglacial meltwater transport. Although the model achieves a good match to empirical evidence without these representations, there is empirical evidence that meltwater transport is an important control on ice velocity evolution (Zwally et al., 2002). Meltwater transport may prove to be more important in other ice sheets contexts, and may improve the representation of the Celtic Sea Ice Stream. Meltwater transport is an important process that remains either not included or highly idealised in ice sheet models.
This work demonstrates that the palaeo record is an ideal test for new ice sheet models that simulate ice stream evolution without the need for statistical inversion. Simulations of palaeo ice sheets have already been compared to RSL records through GIA modelling (Simpson et al., 2009;Kuchar et al., 2012;Auriac et al., 2016), and flowsets (Patton et al., 2016). In addition, we suggest that new ice sheet models developed with the aim of projecting the evolution of the Greenland and Antarctic Ice Sheets over the centennial and millennial timescale should be tested against empirical data of ice stream flow direction from well-constrained sectors of the palaeo record to test for sufficient model skill.

Conclusions
We developed the BISICLES ice sheet model to simulate the spontaneous generation and evolution of ice streams during the advance and retreat of ice sheets. We model basal sliding by using a Coulomb law dependent on basal water pressure, with a simple scheme for local calculating the local balance of water in a till layer. This development allows for the simulation of ice streams that evolve along with the ice sheet. We tested the skill of an ice sheet model's ability to simulate ice streams using the palaeo-record of the BIIS. The palaeo record offers the opportunity to test the skill of the model at simulating ice stream location, spacing, and evolution across a multi-millennial glacial cycle. Despite several assumptions, namely a uniform bed friction coefficient, a highly idealised surface mass balance forcing, and uniform maximum till-water layer thickness, the simulated ice streams strongly match empirically reconstructed ice streams of the BIIS. Results of our simulations suggest that the major controls on ice streaming of the BIIS are topographic troughs and marine margins. We found a spatially variable geothermal heat flux to have only a minor local impact of the pattern of ice streaming, where topographic features are less pronounced.
The simulated ice streams of the BIIS demonstrates the high model skill of ice stream placement and evolution over millennial timescales. With application of this model to more diverse palaeo and modern settings, it could be determined if the model has sufficient skill to model the future evolution of the Antarctic and Greenland ice sheets over the medium-and long-term. Using a model that adequately accounts for ice stream evolution is crucial to future projections of ice sheet change. Ice sheet simulations working in tandem with empirical data therefore can help progress our simulations of future ice sheet evolution towards greater accuracy.

Code and data availability
We used a branch of the BISICLES ice sheet model, revision 3776 (https://anag-repo.lbl.gov/svn/BISICLES/public/branches/slc_dev_ 2018). Output files and required input files to reproduce the described experiments are accessible by contacting the corresponding author and will be deposited independently for final publication.

Author contributions
NG designed the study with LJG and JCE. SLC and LJG developed the model code. NG ran the simulations, performed the analysis, and wrote the manuscript with contributions from all authors. JCE, CDC, and DMH contributed to the comparison to empirical data and interpretation of the results. All authors contributed to the writing of the manuscript.

Conflicts of interest
The authors declare that they have no conflict of interest.