Collapse of the Last Eurasian Ice Sheet in the North Sea Modulated by Combined Processes of Ice Flow, Surface Melt, and Marine Ice Sheet Instabilities

The record of the confluence and collapse of the British‐Irish Ice Sheet and the Fennoscandian Ice Sheet is obscured by the North Sea, hindering reconstructions of the glacial dynamics of this sector of the Eurasian Ice Sheet complex during the last glacial cycle. Previous numerical simulations of the deglaciation of the North Sea have also struggled to capture the confluence and separation of the British‐Irish and Fennoscandian Ice Sheets. We ran an ensemble of 70 experiments simulating the deglaciation of the North Sea between 23 and 18 ka BP using the BISICLES ice sheet model. A novel suite of quantitative model‐data comparison tools was used to identify plausible simulations of deglaciation that match empirical data for ice flow, margin position, and retreat ages, allowing comparisons to large amounts of empirical data. In ensemble members that best match the empirical data, the North Sea deglaciates through the collapse of the marine‐based Norwegian Channel Ice Stream, unzipping the confluence between the British‐Irish Ice Sheet and the Fennoscandian Ice Sheet. Thinning of the Norwegian Channel Ice Stream causes surface temperature feedbacks, rapid grounding line retreat, and ice stream acceleration, further driving separation of the British‐Irish and the Fennoscandian Ice Sheets. These simulations of the North Sea deglaciation conform with the majority of empirical evidence, and therefore provide physically plausible insights that are consistent with reconstructions based on the empirical evidence, and permit a quantitative comparison between model and data.

Despite previous work, the deglaciation dynamics of the North Sea sector remain poorly understood. Indeed, Hughes et al. (2016) suggests that the North Sea and the eastern Barents-Kara Sea, are the "most ambiguous sectors" of the Eurasian Ice Sheet complex. The marine nature of the catchment makes empirical work challenging and studies of the controls of retreat are limited (e.g., Cotterill et al., 2017;Dove et al., 2017;Emery et al., 2019;Phillips et al., 2018;Roberts et al., 2018). With limited empirical data, preliminary reconstructions of the BIIS suggested the ice sheet had extensive marine sections and confluence with the Fennoscandian Ice Sheet (Geikie, 1863). As evidence grew in the second half of the 20th century, reconstructions of the BIIS and Fennoscandian Ice Sheet did not suggest ice sheet confluence in the North Sea (Boulton et al., 1985). However, more extensive offshore work in the 1990s and 2000s revealed glacial sediments in shallow marine cores, and interpretations from seismic reflection surveys across the North Sea (Graham et al., 2007;Sejrup et al., 1994), once again suggested confluence between the two ice sheets. The southern limit of the ice sheet is constrained with a combination of terrestrial and marine data, limits in Norfolk , the Bolders Bank , and Denmark (Houmark-Nielsen & Kjaer, 2003). These limits are consistent with a marine margin across Dogger Bank (Sejrup et al., 2016). Improved bathymetric data has also been used to map ice sheet flow and margin features in the North Sea (Bradwell et al., 2008;Sejrup et al., 2016).
An extensive confluence between the BIIS and Fennoscandian Ice Sheet is now broadly accepted (Graham et al., 2007;Sejrup et al., 1994), but the mechanism and timing of the subsequent separation of the ice sheets is unclear. For example, using evidence from high resolution geomorphological mapping Sejrup et al. (2016) suggested the deglaciation of the North Sea was triggered by the retreat of the large Norwegian Channel Ice Stream (NCIS), leading to the debuttressing of adjacent ice, and an "unzipping" of the BIIS and Fennoscandian Ice Sheet along the Norwegian Channel. However, other studies (Bradwell et al., 2008;Carr et al., 2006) suggest the initial separation of the two ice sheets occurred to the west of the Norwegian Channel, initially retreating back into the deeper bathymetry of the Witch Ground Basin, forming a calving GANDY ET AL.  bay ( Figure 1). In this scenario, the NCIS remained fully advanced to the continental shelf edge while parts of the North Sea deglaciated.
A reconstruction of deglaciation based on the DATED-1 database followed the reconstructed Bradwell et al. (2008) calving bay, suggesting the initial ice sheet separation was to the west of the Norwegian Channel (Hughes et al., 2016). In this reconstruction, a large embayment in the north of the North Sea sector developed, while confluence between the BIIS and Fennoscandian Ice Sheets remained in the south of the North Sea, and the NCIS remained extended to the trough mouth. However, the sparsity of dates in the North Sea means that these reconstructions are tentative (Hughes et al., 2016). For example, there is uncertainty and conflicting evidence of the timing of retreat of the NCIS (King et al., 1998;Sejrup et al., 1994;Svendsen et al., 2015). Additionally, the formation of a large embayment in the north of the North Sea while the NCIS remained extended and active is enigmatic, and potentially glaciologically implausible (Clark et al., 2012;Hughes et al., 2016), as there is no contemporary analog for an ice stream laterally unbounded along one flank. Overall, although empirical constraints on the deglaciation of the North Sea have improved, considerable uncertainty remains regarding both the location of initial separation of the BIIS and Fennoscandian Ice Sheets (Bradwell et al., 2008;Sejrup et al., 2016), and subsequent evolution and speed of ice retreat in the southern North Sea Roberts et al., 2018).
Numerical modeling has not resolved the uncertain deglaciation dynamics of the North Sea either. Boulton and Hagdorn (2006) simulated the confluence of the BIIS and the Fennoscandian Ice Sheet, and found that due to a diminished ablation area after initial confluence the ice sheet thickness increased rapidly to form an ice divide running east-west across the North Sea. This has implications for ice flow, meaning it is unlikely the NCIS flowed the ∼700 km length of the Norwegian Channel concurrently. This numerical modeling work helped to address some uncertainties in regard to the confluence of the BIIS and Fennoscandian Ice Sheets, but does not simulate the subsequent deglaciation of the North Sea sector. Furthermore, these shallow ice approximation simulations did not account for the longitudinal stresses needed to accurately simulate ice stream evolution (Hindmarsh, 2009) or grounding line migration (Schoof, 2007), likely to have been important to the deglaciation of the North Sea.
More recently, a series of simulations investigated the advance (Patton et al., 2016) and retreat (Patton et al., 2017) of the Eurasian Ice Sheet complex. Despite broadly matching empirical data across other sectors of the Eurasian Ice Sheet complex, the simulations are unable to simulate North Sea deglaciation if the maximum extent is consistent with empirical evidence of full glaciation. Therefore, there have been no previous simulations of the North Sea deglaciation that are broadly consistent with empirical evidence. We hypothesize that this is because of the limitations in the climate forcing or ice flow physics previously used.
Recent advances in ice sheet modeling permit new simulations of the North Sea deglaciation. In particular, the L1L2 physics and adaptive mesh of the BISICLES ice sheet model allow accurate simulation of marine sectors of ice sheets (Cornford et al., 2013), such as the northern North Sea. A new basal sliding scheme coupled with a hydrology parameterization also allows the simulation of spontaneous ice stream generation and evolution during ice advance and retreat, with good match to empirical evidence (Gandy et al., 2019). This means that two important controls on the North Sea deglaciation-marine dynamics and the influence of the NCIS-can be accurately simulated for the first time.
Here, we present an ensemble of simulations of the deglaciation of the North Sea, using the state of the art BISICLES ice sheet model. We aim to simulate deglaciation in a manner that conforms with the majority of empirical evidence, and through this determine the likely mechanisms and style of deglaciation of the North Sea. Two research questions, motivated by the considerable uncertainty of deglaciation style in the North Sea, frame the analysis of the work; 1. What was the role of the Norwegian Channel Ice Stream in the deglaciation of the North Sea? 2. How do ice stream dynamics interact with other mechanisms driving deglaciation of the North Sea?

Methods
We use the BISICLES marine ice sheet model to simulate the North Sea sector of the Eurasian Ice Sheet complex. BISICLES is an ice sheet model with L1L2 physics simplified from the full Stokes flow equations (Schoof & Hindmarsh, 2010), using an adaptive mesh to allow the horizontal resolution to be increased at the grounding line, at the margin, and at regions of high velocity (Cornford et al., 2013). This means BISICLES can simulate grounding line migration without parameterization, unlike models previously used to simulate the deglaciation of the North Sea (Patton et al., 2017). BISICLES has previously been successfully used to simulate the evolution of contemporary (Cornford et al., 2016;Favier et al., 2014;Gong et al., 2017) and palaeo (Gandy et al., 2018(Gandy et al., , 2019 ice sheets. We set-up BISICLES in a manner similar to Gandy et al. (2019), in that an idealized basal hydrology scheme is coupled to a Coulomb sliding law dependent on ice and basal water pressure, in order to simulate regions of ice streaming. This method has successfully simulated the majority of ice streams of the BIIS during advance and retreat of the ice sheet (Gandy et al., 2019). For the majority of experiments, we set a 16 × 16 km grid refined once in the North Sea sector to produce a horizontal resolution of 8 × 8 km, except in the high-resolution experiments subsequently described, which have a maximum horizontal resolution of 1 × 1 km.
Initially, a 10,000-year spin-up step simulates the growth and confluence of the BIIS and Fennoscandian Ice Sheet from limited ice extent to extensive glaciation. The spin-up process is described in detail in the supplementary information. Of particular note is that a negative Surface Mass Balance (SMB) anomaly is imposed in the Southern North Sea to help simulated margins to resemble the reconstructed margins. Without this forcing, the margin advances considerably further south than is reconstructed from empirical data, meaning the deglaciation simulations would not start from a reasonable extent, and hence would be unlikely to produce a realistic deglaciation pattern. More details on the implementation and motivation of this correction are provided in the supplementary information.

Deglaciation Ensemble
The end of the spin-up step produces an ice sheet in good agreement with ice extent reconstructed for 23 ka BP (Bradwell et al., 2019;Roberts et al., 2018), which is used as the starting point for a 70-member deglaciation ensemble. Each ensemble member is then run for a 5,000-year deglaciation phase, where seven parameters are varied (Table 1). SMB is calculated every 25 years using the PyPDD model (Seguinot, 2013), forced using contemporary surface temperature and precipitation data (Fick & Hijmans, 2017) (Gordon et al., 2000;Pope et al., 2000;Valdes et al., 2017), and are a refinement of those previously reported by Singarayer et al. (2011) with updated boundary conditions including ice mask, ice orography, bathymetry, and land-sea mask (Ivanovic et al., 2016). An east-west precipitation gradient is imposed using the same pattern employed by Patton et al. (2016Patton et al. ( , 2017) that aims to capture the influence of the Eurasian Ice Sheet complex on precipitation patterns. The modeled ice sheet surface is used to correct surface air temperatures every 25 years through the deglaciation, linearly interpolating based on an environmental lapse rate (γ). Using the same GCM results as used to calculate SMB, subshelf melt rate is derived from Sea Surface Temperature (SST) values, using a linear relationship between SST (T) and subshelf melt rate (ssm) (Rignot & Jacobs, 2002); (1) Using the unadjusted SST values from the GCM produces a subshelf melt rate of 63.2 m y −1 adjacent to the Norwegian Channel Ice Stream at 28 ka BP. This is an idealized parameterization to calculate subshelf melt rates, and in reality, the melt rate is dependent on the stratification of vertical evolution of the water column (e.g., Alvarez-Solas et al., 2019;Ezat et al., 2014). Accurately calculating subshelf melt rates is a significant challenge, and this motivates the idealized approach in tandem with using a wide range of adjustment values in the ensemble. To recreate isostatically adjusted bed topography, we modify modern topography using reconstructions from a Glacio-Isostatic Adjustment model, similar to Bradley et al. (2011), with updated ice-loading based on new margin reconstructions  and resulting ice thicknesses from the ICESHEET ice sheet model (Gowan et al., 2016). GEBCO (Becker et al., 2009) provides modern offshore bathymetry, and SRTM (Farr et al., 2007) provides onshore topography. The RSL change is updated every 1,000 years, and topography is linearly interpolated between these points.

Model-Data Comparison
In order to determine which of the 70 ensemble members simulate behavior consistent with the available empirical data, we use new quantitative model-data comparison tools (Ely et al., 2019a(Ely et al., , 2019b. Model-data comparison tools are used to identify ensemble members that match empirical data. The Automated Proximity and Conformity Analysis (APCA) tool Napieralski et al., 2007) is used to compare modeled ice sheet margins and empirically mapped moraines. A subset of 132 moraines larger than 32 km long is used Dove et al., 2017;Sejrup et al., 2016). Smaller moraines are ignored because it is unlikely the coarse model resolution could provide a conformity match in a meaningful manner. A model-data match is defined when both the proximity is below a threshold of 16 km (two model gridboxes), and conformity is below a threshold of 8 km (one model gridbox). These values are chosen as the ability of the model to match margins is limited by horizontal grid resolution.
The Automated Timing Accordance Tool (ATAT) (Ely et al., 2019a(Ely et al., , 2019b) is used to identify matches between modeled deglaciation ages and geochronological data (in calibrated years before present). Geochronological dates with a quality control rating of Green or Amber by Small et al. (2017) were used, along with additional offshore dates from Bradwell et al., 2021;Callard et al., 2018;Evans et al., 2018Evans et al., , 2019Evans et al., , 2021Roberts et al., 2018Roberts et al., , 2019, resulting in 131 geochronological dates for comparison with modeled deglaciation ages. ATAT calculates the wRMSE between modeled deglaciation ages and geochronological data, accounting for the uneven spatial distribution of dates. ATAT also calculates the percentage of ice-free dates agreed within error, which is used to rank the ensemble members.
The Automated Flow Direction Analysis (AFDA) tool  is used to identify matches between mapped flowsets, which are derived mostly from drumlins, and modeled flow direction. Flowsets larger than 256 km 2 (4 gridboxes at model resolution) are taken from Greenwood and Clark (2009) and Hughes et al. (2014), resulting in 83 flowsets for comparison with the modeled flow directions. AFDA calculates the mean residual angle and variance of offset between modeled and empirically derived flow directions. A model-data match is when both the mean residual vector within a flowset is below a threshold of 15°, and the mean variance is below a threshold of 0.06. These values are more lenient than previously set (e.g., Ely et al., 2019aEly et al., , 2019b, due to the coarser horizontal resolution, but matching of modeled and mapped flow directions within 15° still provides a strong constraint on ice geometry and flow. AFDA can also determine how often cross-cutting flowsets are replicated with the correct relative timing. This is a difficult test for ice sheet models, and previous applications to the BIIS did not match any cross-cutting relationships (Ely et al., 2019a(Ely et al., , 2019b. Given this, and the reduced sample of cross-cutting once smaller flowsets are removed, this feature of AFDA was not used here.
We define Not Ruled Out Yet (NROY) simulations as those that match more than 75% of moraines using APCA, more than 60% of dates using ATAT, and 50% of flowsets using AFDA. These thresholds were decided through an analysis of the distributions of scores (Supporting Information S5) and exercising our expert judgment of what data we would not expect the model simulations to match. Matching moraines is the least stringent of the comparisons (Ely et al., 2019a(Ely et al., , 2019b, hence the higher threshold. Conversely, some deglaciation ages and flowsets are difficult for the model to match in any simulations, particularly smaller flowsets, so a more lenient threshold is set.

Ensemble Results
The ice advance starts with separate ice sheets covering the BIIS and Fennoscandian Ice Sheet (Figures 2a  and 2b). After 6,000 years of advance, there is confluence between the BIIS and the Fennoscandian Ice Sheet. The initial confluence occurs north of Dogger Bank, while the NCIS is already well advanced in the Norwegian Channel. The straight southern margin in the North Sea is caused by the SMB anomaly imposed to allow the margin extent to resemble empirical reconstructions.
In the deglaciation ensemble, the majority of ensemble members (59/70) simulate the separation of the BIIS and Fennoscandian Ice Sheets. In those simulations, the final separation of the BIIS and the Fennoscandian Ice Sheet occurs in the southern North Sea, near the Jutland Bank ( Figure 1). All 59 simulations that have ice sheet separation have a common deglaciation pattern. Beginning with full confluence of the BIIS and Fennoscandian Ice Sheet (Figures 3a and 3b), there is an initial retreat of the NCIS, with ice remaining to the east and west of the Norwegian Channel. As the grounding line of the NCIS continues to retreat in the Norwegian Channel, there is limited margin retreat to the west in the North Sea. In all simulations, ice cover remains on the Shetland Islands while the NCIS has considerably retreated (Figures 3e-3h).
In the majority of ensemble members, considerable NCIS retreat is concurrent with retreat of the western portion of the BIIS, resulting in limited ice cover in Ireland and a retreat of the Barra Fan Ice Stream (Figures 3g-3j). Meanwhile, margin retreat in the southern North Sea and the Fennoscandian Ice Sheet is limited. This deglaciation pattern results in a narrow marine-based ice cap forming north of the Dogger Bank (Figures 3k and 3l). This narrow cap is slow flowing and entirely recedes within 1,000 years.
The point of separation between the BIIS and the Fennoscandian Ice Sheet is near the Jutland Bank (Figure 1) and is consistent across the 59 ensemble members in which the separation occurs. The legacy of this separation style is that the modeled BIIS ice is often extensive in the North Sea, while significant deglaciation of the western margin of the BIIS has already occurred. This is consistent with the deglaciation pattern reconstructed along the Yorkshire and Lincolnshire coast (Bateman et al., 2015;Clark et al., 2012), and the required damming of proglacial lakes like Glacial Lake Pickering , explained through a North Sea Lobe feature (Bateman et al., 2015;Sutherland et al., 2020). The deglaciation style in the majority of ensemble members causes a Dogger Bank ice dome (Figures 3g-3l). Ice streaming down the Yorkshire and Lincolnshire coast occurs in the majority (64/70) of simulations, and in 22 simulations this ice stream terminates with a lobate margin.
A simulated Dogger Bank ice dome still covers a large portion of the North Sea in the later stages of deglaciation (Figures 3g-3l), and therefore ice remained north of the Dogger Bank for 1-2 ka. By the end of the deglaciation simulation (at 18ka BP), the North Sea is fully deglaciated in 11.9% of the simulations, there is no British Isles ice cover in 7% of simulations, a small Scottish Ice Cap in 12.6% of simulations, and Fennoscandian ice cover in all simulations.

Model-Data Comparison
As expected, the quantitative model-data comparison tools show considerable variations in scores of model-data match (Figure 4)  a best matching ensemble member matching 94.6% of moraines. Of the entire ensemble, 39 simulations matched more than 75% of moraines.
ATAT provides a value of both the percentage of dates where model-data agreement occurs, and the wRM-SE of model-data difference for ice-covered dates where model-data agreement occurs. The wRMSE values are important because it is possible to "match" a high percentage of deglaciation ages, but have a high error on each age. The wRMSE values calculated here are low compared to uncertainties in the empirical data, with a mean wRMSE of 816 years (and a maximum of 1,023 years). In comparison, the errors on a Terrestrial Cosmogenic Nuclide date is in the order of 800 years (Small et al., 2017). On average, each ensemble member matches 52.7% of dates, with the poorest matching ensemble member matching 16.7% of dates, and the best matching ensemble member matching 74.5% of dates. The poorest matching ensemble members do not deglaciate sufficiently, resulting in the majority of geochronological dates remaining ice covered at the end of the simulation. Some early deglaciation ages in central Scotland are never matched while also matching the majority of later (post-18 ka BP) deglaciation ages (Figures 4c and 4d). Of the entire ensemble, 20 simulations matched more than 60% of the deglaciation ages.
Finally, the AFDA tool provides a percentage match between empirically mapped flowsets and modeled flow directions. This is the most stringent of the three quantitative tests because flow direction depends on a number of uncertain controls, such as the ice sheet thermal structure, the ice margins, evolution history, and the palaeo topography. The average match over the ensemble is 43% of mapped flowsets, with the poorest ensemble member matching 16% of flowsets, and the best scoring ensemble member matching 70% of flowsets. Of the entire ensemble, 24 simulations matched more than 50% of flowsets. Of these, 17 simulations matched more than 75% of margins and 60% of deglaciation ages too, meaning that these 17 simulations form the NROY subset of simulations ( Figure 5).
There is no clear distinction between the parameter values of the NROY simulations and those of the rest of the ensemble ( Figure S4.1), indeed the range of values for each of the seven parameters is almost the same in the full ensemble as within the subset of NROY simulations. This is partly because all parameters are varied in tandem and the parameter effects can compensate each other. As an example, a simulation with low precipitation (causing rapid retreat) may also have a high lapse rate (slowing retreat), meaning that NROY simulations appear to be unconstrained in the parameter space. In particular, the full range of values for PDD corrections, lapse rates, Weertman friction coefficients, and Till water depths are represented in the NROY ensemble members. Some limited regions of the parameter space are never or rarely represented in the NROY simulations, such as the lowest values of the Coulomb friction angel (i.e., <0.44), high subshelf melt rates, and negative precipitation corrections. Gaussian process emulation of the model output as a function of the input parameters could be used to investigate the effect of individual parameters on ice sheet behavior.
GANDY ET AL.

High-Resolution Snapshot Results
The initial 70 ensemble members have a minimum horizontal resolution of 8 km and provide evidence of rapid retreat periods in the North Sea sector. However, for grounding line dynamics to be accurately simulated, the horizontal resolution at the grounding line must be 1 km or finer (Gladstone et al., 2012). One ensemble member from the NROY simulations (ns_016) has two periods of rapid grounding line retreat of the NCIS, behavior which is also represented in the other NROY simulations. We re-simulate these two periods of rapid retreat using 1 km horizontal resolution at the grounding line to better capture the dynamics of the rapid retreat. The BISICLES mesh refinement allows 1 km horizontal resolution within the regions shown in Figures 6b and 6c, and not for the rest of the domain. Other than the higher horizontal resolution, the experimental design is identical to the design of the coarser resolution ensemble members.
The simulations are rerun from 21.9 and 20.6 ka BP, for 320 years and 240 years respectively. The model mesh is refined down to 1 km at the grounding line at the regions shown in Figures 6b and 6c. The results GANDY ET AL.
10.1029/2020JF005755 11 of 18 show rapid grounding line retreat in regions of retrograde bed slope (Figures 6d and 6e), indicative of Marine Ice Sheet Instability. The simulated retreat is rapid; in both snapshot simulations, the grounding line has a maximum retreat rate of 2-3 m y −1 (Figures 6d and 6e). This grounding line retreat rate is comparable to simulations of future Thwaites Glacier grounding line retreat (Cornford et al., 2015). Both simulations, and the later simulation in particular, show considerable thinning of the ice sheet upstream of the grounding line during the period of rapid retreat. In the later retreat period, there is upstream thinning of ∼2-3 m y −1 through the rapid retreat, which is comparable to observed thinning rates of the contemporary Pine Island Glacier (Shepherd et al., 2001;Wingham et al., 2009).
The behavior of the high-resolution simulations and the lower resolution ensemble member is broadly consistent, in that it demonstrates a period of rapid retreat, starting and concluding in the same regions of the Norwegian Channel. The simulated grounding line retreats at the same rate, though in most areas the grounding line retreats slightly further (<10 km) in the finer resolution simulations ( Figure S6). This provides confidence in the validity of the periods of rapid retreat simulated in the other ensemble members. However, using high-resolution simulations, more detailed behavior of the grounding line is simulated.

Combined Instabilities
The simulations show the fast retreat of the NCIS is driven by the interplay of ice flow, surface lowering, and marine ice sheet instabilities. At the beginning of the deglaciation ensemble, with full ice extent over the North Sea, the NCIS has an effect on the surface profile of the ice sheet, causing a large (∼500 m) depression in the surface of the ice sheet (Figure 7a). This is similar to the surface expressions of similar contemporary ice streams, like Pine Island Glacier and Thwaites Glacier (Howat et al., 2019).
This large surface depression is important because it causes feedbacks between ice dynamics, further surface lowering, and ice geometry. The initial depression means that at the start of the deglaciation experiments the NCIS was vulnerable to surface melting (Figure 7b), allowing the surface of the NCIS to be lowered further. This has four effects. First, the lowering of the ice stream surface compared to neighboring ice allows ice stream acceleration (Robel & Tziperman, 2016), driving further lowering. Second, the thinning exacerbates the lapse rate feedback that allows further surface melting as the surface lowers. Third, any thinning of the ice stream (which is located in a deep marine trough) means the ice stream becomes more vulnerable to rapid retreat as it approaches buoyancy. Finally, any grounding line retreat into a region of retrograde bed slope can cause marine ice sheet instability. The NROY simulations spanning a wide area of the parameter space ( Figure S4.1) suggests these interacting mechanisms are the primary influence on deglaciation.
Notably, the NCIS remained grounded past the Skagerrak Straight by 18 ka BP in only four ensemble members. Those simulations show an advance of the BIIS and Fennoscandian Ice Sheets in other sectors. In the remaining ensemble members, the combined instabilities of ice flow, surface elevation feedbacks, and marine influence mean that once the NCIS begins to retreat it is likely to continue to retreat fully. However, the NCIS does not experience an uninterrupted retreat for the entire length of the Norwegian Channel. There are periods of rapid grounding line retreat, and periods of grounding line stability. This behavior is supported by the mapping of grounding zone wedges in the Norwegian Channel (Morén et al., 2018) and by the variable ice rafted debris fluxes along the northern North Sea margin (Becker et al., 2018).
This variability of ice flux is evident in the simulations even with a comparatively smooth climate and ocean forcing. Including a more temporally variable atmosphere and ocean forcing, as used in previous simulations (e.g., Hubbard et al., 2009;Patton et al., 2017), would further increase the variability of retreat rate. Warm water intrusion into the Norwegian Channel would be sensitive to relatively small changes in the vertical thermal structure of the ocean and relative sea level. The sensitivity of the Bjørnøyrenna Ice Stream (a large marine ice stream between Norway and Svalbard) to subshelf melt forcing has been demonstrated (Petrini et al., 2018), and given the relative proximity of the NCIS and Bjørnøyrenna Ice Stream it is likely both ice streams were subject to similar variability of marine forcing.

Styles of Deglaciation
As discussed, empirical reconstructions of the separation of the BIIS and Fennoscandian Ice Sheet do not provide a consensus on the style of deglaciation. Bradwell et al. (2008) suggested that the separation began with the formation of a calving bay to the west of the Norwegian Channel, forming a large embayment in the northern North Sea. In contrast, Sejrup et al. (2016) suggested that deglaciation of the North Sea begun with the initial retreat of the NCIS, which debutressed the BIIS to the west, causing no initial significant margin retreat of the BIIS in the northwest while the NCIS retreated significantly, leading to an "unzipping" between the BIIS and the Fennoscandian Ice Sheet. The ensemble of simulations presented in this manuscript allows us to test the plausibility of each deglaciation scenario.
Here, none of the ensemble members began deglaciation of the North Sea through the formation of an embayment in the Witch Ground Basin while the NCIS remained extended. Instead, there are key similarities between the "unzipping" scenario and the majority of simulations described here. In all simulations, the retreat of the NCIS precedes retreat into the Witch Ground Basin. In the 17 NROY simulations, the NCIS retreats past the Skagerrak Straight ( Figure 1) before BIIS North Sea Ice has entirely receded ( Figure 5). Although the Witch Ground Basin is a major topographic depression in the context of the relatively flat North Sea, it is a minor depression in comparison to the neighboring Norwegian Channel, and no simulation in our ensemble demonstrates ice retreating in the Witch Ground Basin but remaining stable in the Norwegian Channel. Previously, it has been suggested that because the NCIS was more vulnerable to marine and GANDY ET AL.

10.1029/2020JF005755
13 of 18 climate forced retreat, initial retreat into the Witch Ground Basin while the NCIS remained fully extended is "enigmatic" (Clark et al., 2012;Hughes et al., 2016). Our findings also suggest that this deglaciation style is unlikely.
In contrast, there are similarities between the style of deglaciation in the NROY simulations and the deglaciation reconstructed by Sejrup et al. (2016). Both the reconstruction and the simulations presented here suggest a key role of the NCIS, initially retreating before the rest of the North Sea. However, while simulations here place the final point of confluence before separation as in the Jutland Bank region, Sejrup et al. (2016) place it further North, to the west of the Norwegian Channel.
The simulations presented here were computationally costly, therefore only the deglaciation was included in the 70-member ensemble, rather than a full glacial cycle. Starting these ensemble experiments from the end of a single advance and confluence simulation rather than a stable maximum extent in part helps represent the transient nature of the last glacial cycle, but it would still be preferable to run an ensemble of transient advance and retreat simulations, if computational cost permitted.
With the current experiments, it is not clear how sensitive the deglaciation style is to different styles of ice sheet build-up. As computational costs reduce, an ensemble of a transient advance and retreat experiments would be advantageous. Alternatively, a computationally cheaper model can be used to simulate a full glacial cycle of the Eurasian Ice Sheet complex (Patton et al., 2016(Patton et al., , 2017, but so far such models are not able to simulate deglaciation of the North Sea in a manner consistent with the empirical record.

Western North Sea Deglaciation
A combination of geomorphological and geological evidence, and dating has resulted in the reconstruction of a so-called North Sea Lobe, flowing south from the Firth of Forth in Scotland and nearly parallel to the east English coast and extending down to the north Norfolk coast Evans & Thomson, 2010;Sutherland et al., 2020). Specifically, the presence of such a North Sea Lobe helps to explain dates showing ice free conditions across central England, while deglaciation is later along the Yorkshire and Lincolnshire coast (Bateman et al., 2015). Such an ice margin is also necessary to dam reconstructed proglacial lakes in northern and central England (Bateman et al., 2015;Davies et al., 2019;Evans & Thomson, 2010).
Although ice streaming down the Yorkshire and Lincolnshire coast is simulated in numerous ensemble members, a feature as prominent as the empirically reconstructed North Sea Lobe is never simulated. Instead, in the majority of simulations a Dogger Bank ice dome forms during deglaciation, where ice remains in the North Sea, and retreating roughly radially from the point where the BIIS and the Fennoscandian Ice Sheets disconnect. This dome is the remnants of the ice divide that forms across the North Sea at maximum extent, and it is what remains when the majority of surrounding ice has retreated. Compared to empirical reconstructions of the last stages of deglaciation in the North Sea (Figures 8a and 8b), the simulated ice margin extends further east toward the Norwegian Channel (Figure 8c Clark et al. (2012), split into (a) Scenario one with a southward flowing North Sea lobe down the east coast and which retreated to the north, and (b) Scenario 2 with retreat mostly directed toward the west following separation of the two ice sheets and loss of a southern North Sea ice dome, and (c) compared to the simulated Dogger Bank ice dome in our ice sheet modeling.
A North Sea ice dome would require a low equilibrium line and should be evident from relative sea level records around the North Sea (Hughes et al., 2016). Topographically or climatically, the southern North Sea is not an intuitive area to reconstruct an ice dome. Despite this, the empirical data in the North Sea is sufficiently sparse to not exclude these simulations from the NROY subset. The simulations are likely not reproducing the actual deglaciation pattern, but cannot entirely be ruled out by the data.

Conclusions
We completed an ensemble of simulations of the deglaciation of the North Sea sector of the Eurasian Ice Sheet complex, using the BISICLES ice sheet model, simulating deglaciation in a manner that conforms with the majority of empirical evidence for the first time. We used a suite a quantitative model-data comparison tools to compare simulations to a large amount of empirical evidence. The simulations show that the Norwegian Channel Ice Stream was influential in the separation of the British-Irish Ice Sheet and the Fennoscandian Ice Sheet, retreating through the interaction of Marine Ice Sheet Instability, elevation-lapse rate effects, and ice dynamics feedbacks. The retreat of the Norwegian Channel Ice Stream has effects on the surrounding ice masses, and facilitates further retreat of the BIIS. The simulations are consistent with a reconstructed style of deglaciation beginning with retreat of the NCIS, and suggest retreat to the west of the Norwegian Channel while the NCIS remains extended is unlikely.
In the later stages of deglaciation, all NROY simulations produce a Dogger Bank ice dome, where ice remains extensive in the North Sea while the BIIS has significantly retreated in the west. We find this to be a plausible alternative to the North Sea lobe, and hypothesize that a remnant Dogger Bank ice dome formed as a consequence of deglaciating an ice divide in the North Sea. The progress made here in simulating deglaciation dynamics that conforms with the majority of empirical evidence is possible because of improved simulation of marine ice sheet dynamics and the evolution of ice streams. This demonstrates the importance of using an ice sheet model of sufficient skill when simulating the future evolution of marine based ice sheets.