New high‐resolution multibeam data in the Gulf of Bothnia reveal for the first time the subglacial environment of a Bothnian Sea Ice Stream. The geomorphological record suggests that increased meltwater production may have been important in driving rapid retreat of Bothnian Sea Ice during deglaciation. Here we apply a well‐established, one‐dimensional flow line model to simulate ice flow through the Gulf of Bothnia and investigate controls on retreat of the ice stream during the post‐Younger Dryas deglaciation of the Fennoscandian Ice Sheet. The relative influence of atmospheric and marine forcings are investigated, with the modeled ice stream exhibiting much greater sensitivity to surface melting, implemented through surface mass balance and hydrofracture‐induced calving, than to submarine melting or relative sea level change. Such sensitivity is supported by the presence of extensive meltwater features in the geomorphological record. The modeled ice stream does not demonstrate significant sensitivity to changes in prescribed ice stream width or overall bed slope, but local variations in basal topography and ice stream width result in nonlinear retreat of the grounding line, notably demonstrating points of short‐lived retreat slowdown on reverse bed slopes. Retreat of the ice stream was most likely governed by increased ice surface meltwater production, with the modeled retreat rate less sensitive to marine forcings despite the marine setting.

1 Introduction Marine‐terminating glaciers are a major source of mass loss from the contemporary Greenland and Antarctic Ice Sheets [Van den Broeke et al., 2009; Cook et al., 2014; Murray et al., 2015], sensitive to both atmospheric and marine forcings, and thus particularly susceptible to climatic and environmental changes. Submarine melting is often considered one of the most important drivers of marine‐terminating glacier retreat due to the potential to drive ice shelf debuttressing and consequent dynamic thinning [Joughin and Alley, 2011; Luckman et al., 2015]. Where the bed increases in depth upstream this can lead to a runaway effect referred to as marine ice sheet instability [Weertman, 1974; Schoof, 2007], with an increase in grounding line depth leading to increased ice flux due to thicker ice and increased calving. Changes in both sea level and tides can influence the depth of the grounding line, with tides further acting to modulate ice velocities and hydrostatic backstress [Anandakrishnan et al., 2003; Thomas, 2007]. Local topography can also control grounding line retreat, with reverse bed slopes increasing the risk of instability, while topographic highs can create pinning points for stabilization. In addition, the sides of troughs can exert a lateral drag, reducing ice flow and slowing retreat [Whillans and van der Veen, 1997], with some evidence for a narrowing in ice stream width aiding short‐lived slowdown of retreat on reverse bed slopes [Gudmundsson et al., 2012; Jamieson et al., 2012]. In addition to marine forcings, the stability of marine‐terminating glaciers can also be strongly influenced by atmospheric warming [Lea et al., 2014a; DeConto and Pollard, 2016]. Changes in accumulation and ablation can induce thinning promoting a dynamic response, while meltwater reaching the ice‐bed interface can enhance basal lubrication and alter the efficiency and topology of the subglacial hydrological system, further perturbing ice dynamics in at least the subannual scale [Zwally et al., 2002; Sole et al., 2013; Doyle et al., 2014]. Surface meltwater can also induce hydrofracture of crevasses leading to increased rates of calving [Benn et al., 2007; Colgan et al., 2016] and can indirectly cause increased submarine melting through the release of supraglacially originating plumes of subglacial meltwater at glacier termini [Motyka et al., 2003; Chauché et al., 2014]. The relative contribution of these individual external forcings on the dynamics and stability of marine‐terminating glaciers remains poorly constrained [Nick et al., 2009, 2010, 2012; Murray et al., 2010; Straneo et al., 2013], providing major uncertainty in understanding their stability and retreat, compounded by the interlinked nature of and nonlinear feedback between forcings. To better predict future response of these catchments, it is imperative to better understand the physical processes and local factors determining grounding line stability and retreat behavior in a diverse range of settings. The record of well‐constrained observations of contemporary marine‐terminating catchments is short [Kjeldsen et al., 2015; Murray et al., 2015] and only informs us of ice flow behavior on the order of hours to decades [Joughin et al., 2010]. However, the geological record of these catchments’ dynamics can improve our understanding of the roles of climate, marine processes, and topography on marine ice stream retreat, particularly over centennial to millennial time scales [Winsborrow et al., 2010; Jakobsson et al., 2011; Jamieson et al., 2014; Jones et al., 2015]. In particular, landform assemblages created by past ice flow regimes and grounding line processes provide important constraints on physical models of glaciological processes. This work presents a combined geomorphological and modeling approach to examining the controls on the retreat and dynamics of a paleo‐ice stream in the Gulf of Bothnia. We investigate this setting because (i) the epicontinental, approximately 100,000 km2 basin contrasts with the ocean‐facing, continental shelf trough environments that have hitherto been the focus of marine ice instability investigation; (ii) the retreat pattern, controls, and effect of this major ice flow corridor on southern Fennoscandian Ice Sheet (FIS) deglaciation are currently poorly understood; and (iii) the availability of excellent new geomorphological data constrains its behavior. Here we apply a well‐established numerical flow line model [Vieli and Payne, 2005; Nick et al., 2010] to investigate the relative sensitivity of Bothnian Sea Ice Stream retreat to both atmospheric and marine forcings. Specifically, we examine the effects of increased ice surface melting, submarine melting at the terminus, and sea level change, in addition to sensitivity to basal topography and catchment geometry. A newly reported, rich landform record of deglaciation [Greenwood et al., 2016], coupled with a physical exploration of ice sheet retreat processes, provide an opportunity to gain an improved understanding of the balance of processes governing a major marine ice sheet catchment and ultimately its stability.

2 The Geomorphological Record of Paleo‐Ice Flow The Gulf of Bothnia is a shallow, epicontinental basin that lies in the heart of the terrain formerly covered by the FIS (Figure 1). The Bothnian Bay subbasin, in the north, is separated from the larger Bothnian Sea to the south by a shallow sill; the Bothnian Sea is similarly separated from the Baltic Sea by the Åland archipelago. These basins are thought to have hosted a variety of glaciodynamic environments throughout the evolution of the FIS: an interior position close to the main ice divide [Kleman et al., 1997], the head of an extensive Baltic Ice Stream [Holmlund and Fastook, 1993; Boulton et al., 2001], and a major corridor of marine‐ and lacustrine‐based deglaciation [De Geer, 1940; Strömberg, 1989; Lundqvist, 2007]. From the Younger Dryas and subsequent deglaciation, a series of ice lobes are conventionally considered to have crossed the Gulf of Bothnia southeastward into Finland, well defined by glacial lineations, interlobate esker corridors, and large lobate terminal moraines [Punkari, 1980; Johansson et al., 2011]. This whole sector deglaciated under marine (periodically lacustrine) conditions [Björk, 1995; Andrén et al., 2011]. The marine limit on the Swedish High Coast is 286 m above present‐day sea level [Berglund, 2004]. Figure 1 Open in figure viewer PowerPoint Baltic Sea Hydrographic Commission, 2013 2014 Stroeven et al. [ 2016 Hughes et al., 2016 (a) Bothnian Sea ice flow line, based on geomorphological interpretations and estimated maximum and minimum ice stream widths (solid, dashed, and dotted white lines, respectively). The black tickmarks at 100 km intervals are placed along the flow line; the green portion of the line denotes the segment for which there is geomorphological evidence for ice streaming velocities; the yellow rectangles mark the identified sticky zones; and the red stars mark the sites of shoreline displacement data from (north to south) Norrbotten, Ångermanland (Swedish High Coast), Gästrikland, and Södertörn. Areas for which multibeam bathymetry data were available are outlined in grey. The displayed background bathymetry is from the Baltic Sea Bathymetry Database [] and topographic compilation from viewfinderpanoramas.org []. The varve‐based deglacial chronology from. [] between 11.6 and 10.1 ka B.P. is depicted in black on the left, with the hypothesized Younger Dryas (11.6 ka B.P.) limit extended across the Baltic Sea. The inset shows the Baltic‐Bothnian Basins in the context of the Fennoscandian Ice Sheet (including the 21 ka B.P. areal ice extent [after.,]). Extracts from the multibeam data illustrating mega‐scale glacial lineations, indicative of the ice stream pathway flowing SSW in the upper part of the (b) flow line and (c and d) SSE downstream. Large and well‐connected meltwater channel systems overprint and incise the ice stream assemblage (Figures 1b and 1c). All reconstructions of paleo‐ice dynamics in the Bothnian Basin have hitherto been drawn from evidence from the present‐day terrestrial domain; direct glacial geological or geomorphological evidence from the Gulf itself has been almost entirely lacking, and the dynamics of ice retreat are thus poorly understood. New, extensive, high‐resolution multibeam echo‐sounding data allow direct investigation of its glacial geomorphological record for the first time, revealing a paleo‐ice stream directed SSW to southward through the Bothnian Sea [Greenwood et al., 2015, 2016]. Mega‐scale glacial lineations (MSGLs) >20 km in length define a slightly sinuous, ~350 km long ice stream pathway, with an onset zone over the Västerbotten coast of Sweden marked by distinct transitions in lineation size and elongation [Greenwood et al., 2015]. The width of this onset zone is limited to ~40 km, while the ice stream trunk, indicated by increasing downstream MSGL length and elongation, continues to a position in the central‐southern Bothnian Sea. Here MSGLs overprint each other in two to three subparallel sets that together indicate a splayed, lobate terminus. Beyond this position, small (55–700 m long) crag and tails and streamlined crystalline bedrock exhibit convergent flow into the Åland Sea. We interpret initial retreat of moderate‐to‐slow flowing ice across the Åland sill. A late‐stage streaming event occurred as ice retreated through the Gulf of Bothnia, but which did not extend the full length of the basin and likely did not fill its entire breadth [Greenwood et al., 2016]. Terrestrial‐based chronologies for deglaciation provide an available time window for complete deglaciation of the gulf of only approximately 1000–1200 years [Hughes et al., 2016; Stroeven et al., 2016]. This marks a twofold to fourfold increase in the ice sheet retreat rate relative to the deglaciation up until the Younger Dryas. Following retreat from the Younger Dryas position, a notable lack of grounding line deposits suggests that there were no major standstills during deglaciation [Greenwood et al., 2016]. The lobate and internally crosscutting arrangement of distal MSGLs suggests the margin briefly paused at the terminus of the Bothnian Sea Ice Stream, where the flow direction close to the grounding line shifted back and forth. Approximately 160 km up ice, a small group of transverse wedges suggest localized grounding line pauses; otherwise, De Geer moraines only appear in the present‐day coastal and terrestrial domains, following loss of ice from the offshore sector. An extensive meltwater landform record comprising interconnected channels, eskers, and large (kilometer‐scale) erosional corridors suggest that meltwater access to and drainage via the subglacial system was plentiful during deglaciation. These geomorphological observations indicate that retreat of the Bothnian Sea Ice Sheet catchment was associated with episodically high ice flow velocities, rapid margin retreat, high surface melting, and access of supraglacial meltwater to the bed. Previous modeling of the FIS has largely focused on reproducing areal ice extent, and to a lesser extent ice dynamics, at the ice sheet scale [Arnold and Sharp, 2002; Forsström et al., 2003; Siegert and Dowdeswell, 2004; Clason et al., 2014]. A number of modeling attempts have reported difficulty in avoiding overly rapid deglaciation through the Baltic‐Bothnian sector of the FIS [Holmlund and Fastook, 1993; Clason et al., 2014], with Holmlund and Fastook [1993] resorting to the prescription of a “sticky spot” of cold‐based ice over the island of Åland (Figure 1) to decelerate ice sheet thinning and dam ice flow. Here we apply a flow line model specifically to the Bothnian Sea case to investigate the possible controls on the retreat behavior indicated by geomorphological evidence. We define an 835 km long central flow line in accordance with the MSGL assemblage and the small‐scale lineations converging on the Sea of Åland. We specify a late Younger Dryas terminal position (11.6 ka B.P., based on the reconstructed margins of Stroeven et al. [2016]) and extend our flow line headward to the approximate position of a late Younger Dryas ice divide [Kleman et al., 1997]. We define a minimum width based on landform evidence of the paleo‐ice stream lateral margin and allow a plausible maximum width to account for broader, nonstreaming episodes and uncertainties due to restricted data coverage. Within this model setup, the defined flow line path is fixed for the duration of each model simulation. We consider the geometry of the modeled ice stream reasonable for the period of time in which ice is in the present‐day offshore and nearshore domain. However, it is well known that the FIS ice divide migrated west during this time and that the retreating flow line would also have swung westward [Hughes et al., 2016; Stroeven et al., 2016]. The uppermost part of our flow line should therefore be treated only as an accumulation zone for the distal, offshore ice dynamics and not as an ultimately terrestrial retreat flow line.

3 Numerical Modeling 3.1 Model Description Vieli and Payne, 2005 Nick et al. [ 2010 Nick et al., 2012 Carr et al., 2015 Lea et al., 2014a 2014b Jamieson et al. [ 2012 2014 van der Veen and Whillans [ 1996 Fowler [ 2010 (1) x is the distance along the flow line, h is the ice surface elevation, ρ i is the density of ice, and ρ p is proglacial water density which is held constant at 1000 kg m−3 to reflect the largely fresh or brackish water of the Late Weichselian Gulf of Bothnia [Björk, 1995 Andrén et al., 2000 g is the acceleration due to gravity; U is the vertically averaged horizontal ice velocity; H is the ice thickness; D is the depth of the base of the ice below the surface of the proglacial water body; A S is the basal roughness parameter; W is the glacier width; A is the temperature‐dependent rate factor [Glen, 1955 μ is the basal friction parameter (typically 1 but can be modified to simulate meltwater‐enhanced basal sliding); m is the friction exponent; and v is the strain rate‐dependent effective viscosity, defined as (2) n is the exponent in Glen's flow law. Evolution of the ice thickness along the flow line accounts for changes in glacier width and surface mass balance, a, following (3) q is the horizontal flux of ice through a cross section of the flow line following q = HWU. Surface mass balance is implemented through a linear gradient which changes with elevation from a minimum at sea level to a maximum at the ice divide, with additional bounds set to restrict overly positive or negative values. The velocity boundary condition at the calving front is calculated as (4) We apply a well‐established, one‐dimensional flow line model [], described in full in. [], to the Bothnian Sea case. The model has previously been applied to Greenland outlet glaciers in both contemporary [.,.,] and historical [.,] settings and used for the simulation of a paleo‐ice stream in Antarctica by. []. The model applies a nonlinear effective pressure‐sliding law, with driving stress (left‐hand term) balanced by the longitudinal stress gradient (first term), basal drag (second term), and lateral drag (third term) following] and], aswhereis the distance along the flow line,is the ice surface elevation,is the density of ice, andis proglacial water density which is held constant at 1000 kg mto reflect the largely fresh or brackish water of the Late Weichselian Gulf of Bothnia [.,].is the acceleration due to gravity;is the vertically averaged horizontal ice velocity;is the ice thickness;is the depth of the base of the ice below the surface of the proglacial water body;is the basal roughness parameter;is the glacier width;is the temperature‐dependent rate factor [];is the basal friction parameter (typically 1 but can be modified to simulate meltwater‐enhanced basal sliding);is the friction exponent; andis the strain rate‐dependent effective viscosity, defined aswhereis the exponent in Glen's flow law. Evolution of the ice thickness along the flow line accounts for changes in glacier width and surface mass balance,, followingwhereis the horizontal flux of ice through a cross section of the flow line following. Surface mass balance is implemented through a linear gradient which changes with elevation from a minimum at sea level to a maximum at the ice divide, with additional bounds set to restrict overly positive or negative values. The velocity boundary condition at the calving front is calculated as Nick et al., 2010 Nye, 1957 Benn et al., 2007 d s , is determined by the longitudinal stresses and crevasse water level, following Benn et al. [ 2007 (5) ρ w is the density of meltwater; d w is the depth of water in the crevasse; and R xx is the normal resistive stress which relates to the longitudinal stretching rate, , following Glen's flow law [van der Veen, 1999 (6) This condition balances the difference between the hydrostatic pressure of the ice and proglacial water with the stretching rate at the terminus. The model includes a physically based full‐depth calving criterion [.,], which combines the depths of surface and basal crevasses within a field of closely spaced crevasses [.,]. Surface crevasse depth,, is determined by the longitudinal stresses and crevasse water level, following. [], such thatwhereis the density of meltwater;is the depth of water in the crevasse; andis the normal resistive stress which relates to the longitudinal stretching rate,, following Glen's flow law [ d b , is calculated following (7) H ab is the height above buoyancy, following (8) Tuning the water depth in surface crevasses thus allows for exploration of the effect of increased meltwater production on calving rates and grounding line retreat. The height of basal crevasses,, is calculated followingwhereis the height above buoyancy, following The formation of a floating tongue or ice shelf is permitted when the ice thickness is less than the flotation thickness. A simple implementation of submarine melting is applied both at the grounding line and where ice is floating, decreasing the ice thickness by a fixed amount for all grid points between the grounding line and the calving front. A moving model grid, for which the horizontal grid spacing adjusts at each time step to accommodate the new glacier length, is used to allow for continuous tracking of the grounding line position [Vieli and Payne, 2005]. Values for constants and physical parameters used in the model are given in Table 1. Table 1. Values for Physical Parameters and Constants Applied in the Model Description Value Acceleration due to gravity, g 9.8 m s−2 Density of ice, ρ i 917 kg m−3 Density of proglacial water, ρ p 1000 kg m−3 Density of meltwater, ρ w 1000 kg m−3 Glen's flow law exponent, n 3 Glen's flow law coefficient, A 2.9 × 10−17 Pa−3 a−1 (ice temperature −5°C) Friction exponent, m 3 Basal friction parameter, μ 1 (unless otherwise stated) Grid size (variable) Approximately 500 m Model time step 0.005 years 3.2 Model Inputs The ice streambed profile was extracted along the flow line (Figure 1) at 500 m horizontal intervals, with present‐day bathymetric values obtained from the Baltic Sea Bathymetry Database [Baltic Sea Hydrographic Commission, 2013] and topographic values from available 50 m gridded elevation data [Lantmäteriet, 2015]. Given the high local topographic variability of the bed, spline smoothing was applied to avoid including features unrepresentative of the bay‐wide basal topography. The central flow line elevation profile was used in modeling, rather than width‐averaged elevation, due to uncertainty in prescribing the width of the ice stream and clear data constraint on the central pathway. To account for bed isostatic evolution and sea level forcing throughout the model period (11.6–9.65 ka B.P.), we use empirically derived relative sea level (RSL) reconstructions that span the operational length of the Bothnian Sea Ice Stream. Shoreline displacement data from sites in Norrbotten [Lindén et al., 2006], Ångermanland [Berglund, 2004], Gästrikland [Berglund, 2005], and Södertörn [Hedenström and Risberg, 1999] were selected accordingly (Figure 1). The published 14C ages for the isolation of lake basins at each site were converted to calendar years B.P. using the IntCal13 calibration curve [Reimer et al., 2013]. Second‐ or third‐degree polynomials were applied to interpolate between the calibrated dates for each site and generate a local RSL curve with a temporal resolution of 10 years. A spline was then fitted to the four shore displacement values at each time step to generate an evolving RSL model that is appropriate in both time and space at each step along our flow line (Figure 2). Figure 2 Open in figure viewer PowerPoint (a) Grounding line retreat for experiments testing sensitivity to the inclusion and exclusion of RSL change, changes to ice stream width (ISW), and changes to bed tilt (BT). The late Younger Dryas grounding line position is shown to be insensitive to RSL change alone (black dashed line). Retreat for ISW and BT tests is forced by surface mass balance scenario SMB_040, and the dashed blue line represents the running SMB_40 with no RSL change. (b) Changes to bed elevation along with flow line over a 2000 year model run (one profile produced every 50 years for 40 time steps), forced by RSL change. The grey area represents the less well‐constrained upper portion of the flow line, above present‐day sea level, here and in subsequent figures. The spatial distribution of an initial surface mass balance was chosen to reflect the colder and drier conditions of the Younger Dryas [Carlson, 2013; Muschitiello et al., 2015]. We applied a surface mass balance that varies linearly with elevation, ranging from −0.5 m a−1 m−2 at sea level up to a fixed maximum of 0.25 m a−1 m−2 at 1000 m above sea level (asl), with an associated equilibrium line altitude of 665 m asl. This is broadly comparable with reconstructed Younger Dryas gradients for northern Norway described in Rea and Evans [2007]. With this setup, we investigated atmospheric and marine forcings on the catchment via four model mechanisms: relative sea level at the grounding line, increasingly negative surface mass balance (reduction of ice thickness across the surface profile), increasing the depth to which surface meltwater fills crevasses (hydrofracture trigger), and increasing submarine melt (reduction of ice thickness across the floating portion). 3.3 Model Experiments The model was spun up to an initial stable, late Younger Dryas position consistent with the 11.6 ka B.P. marginal position presented in Stroeven et al. [2016]. This steady state configuration, from which all perturbation experiments were initialized, was produced with sea level held constant at spatially interpolated values for 11.6 ka B.P., a constant crevasse water depth of 95 m (47% of the steady state terminus thickness), an average of the maximum and minimum ice stream widths (Figure 1), and no submarine melting. The surface mass balance was initially set at 0.2 m a−1 m−2 across the full length of the profile for 1500 years to allow for the growth and evolution of an ice surface with a realistic slope profile. A surface mass balance varying linearly with elevation, as described above, was then applied, and the model was run until the ice surface had stopped evolving after 6500 years. Thereafter, our model experiments were run for a period of 2 kyr, reflecting the time between the late Younger Dryas and full deglaciation as described by Stroeven et al. [2016]. Four sets of experiments, a total of 75 runs, were conducted as described in Table 2, designed to explore the sensitivity of ice flow behavior and retreat rates to (i) system geometry, (ii) atmospheric forcings, (iii) marine forcings, and (iv) combined forcings. In each case, the initial variable values were held at those used in the spin‐up (as above). Table 2. Description of Experiment Sets Testing Sensitivity to Changes in System Geometry and Sensitivity to Linear and Stepwise Changes in Atmospheric and Marine Forcings Experiment Set System Geometry Crevasse Water Depth (m) Mass Balance at Sea Level (m a−1 m−2) Submarine Melt (m a−1) Bed Profile Ice Stream Width Linear Increase Over 500 Years From 95 Step Increase Every 500 Years From 95 Linear Decrease Over 500 Years From −0.5 Step Decrease Every 500 Years From −0.5 Linear Increase Over 500 Years From 0 Step Increase Every 500 Years From 0 1. Sensitivity to bed topography, ice stream width, and RSL change Normal, 10% tilt, 20% tilt Average, 50%, 75%, 125%, 150% Fixed at 95 ‐ 0.4 ‐ Fixed at 0 ‐ 2. Sensitivity to individual atmospheric forcings Normal Average 5, 10, 15, 20, 25, 30 5, 10, 20, 30 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.65, 0.9 0.1, 0.2, 0.3, 0.4 Fixed at 0 ‐ 3. Sensitivity to individual marine forcings Normal Average Fixed at 95 ‐ Fixed at −0.5 ‐ 2, 5, 10, 15, 20, 25, 50, 100, 250 10, 25, 50, 100 4. Sensitivity to combined forcings Normal Average 5, 10, 15 ‐ 0.1, 0.2, 0.3 ‐ 2, 5, 10 ‐ Each model run is named in accordance with the following nomenclature. The code prefix indicates the type of forcing applied, whereby BT is the bed tilt, ISW is the ice stream width, SMB is the surface mass balance, CWD is the crevasse water depth, SM is the submarine melt, and BF is the basal friction. In the case of model geometry (BT and ISW), the suffix indicates the deviation from the standard geometry, expressed as a percentage, and for BF the suffix indicates the fraction of full basal friction (=1.0). For SMB, CWD, and SM, each in their respective units, the suffix indicates the rate of change of the forcing per 500 years (Table 2). The default run setup for SMB, CWD, and SM is a linearly applied forcing; experiments using step changes include the term “step” in the suffix. Following this scheme, experiment SMB_step040 denotes forcing the model with surface mass balance (holding all other parameters at their default/initial value), applied in a stepwise fashion with a change of −0.4 m a−1 m−2 every 500 years to a final magnitude of −1.7 m a−1 m−2. The first set of experiments was designed to test sensitivity to topography and geometry, including the isostatic adjustment of the topography manifest in our RSL model. An initial experiment forced retreat solely with RSL change, holding other variables constant; an additional test of RSL sensitivity removed the RSL forcing from a standard (surface mass balance forced) model run with default catchment geometry. Thereafter, sensitivity testing to both bed tilting and ice stream width was investigated, using the same surface mass balance forced run (SMB_040) as a control. Theory dictates that bed depth, and thus slope, is a strong control on the rapidity of ice stream retreat and grounding line stability [Weertman, 1974; Schoof, 2007]. While the spatially nonlinear isostatic rebound in Scandinavia [cf. Hedenström and Risberg, 1999; Berglund, 2004] is accounted for in our RSL adjustment, there are uncertainties inherent in radiocarbon dating and calibration of samples used to infer RSL history, and further errors are likely introduced in interpolating between sites. Sensitivity testing to bed tilt was implemented by increasing the slope of our RSL‐adjusted bed topography by 10% and 20%. We explore sensitivity of the model to the defined catchment width by increasing and decreasing the average flow path width (model default) by 25 and 50% to cover the drawn extent of maximum and minimum ice stream boundaries. The influence of atmospheric warming and ice surface melting is investigated via changes in surface mass balance and crevasse water depth in the second set of experiments. For both we test the response to linearly increasing forcings and stepwise changes in forcing. Surface mass balance is applied along the flow line, forcing evolution of the ice surface in combination with ice flux, and changes in ice stream width (equation 3). Crevasse water depth is implemented as d w in equation 5, whereby increasing crevasse water depth allows for deeper propagation of surface crevasses via hydrofracture and the possibility for increased calving. Crevasse water depth is a poorly constrained absolute value and in this implementation should not be treated as a real value but rather a tuning parameter linked to a physical mechanism. An initial water depth of 95 m produces a stable spin‐up state, so we begin increasing values from here. Linear forcing experiments impose a steady yearly increase, while step change experiments apply the specified forcing once every 500 years (Table 2). For CWD, values increase from an initial 95 m to an end‐of‐run value ranging from 115 (CWD_5) to 215 m (CWD_30). In the case of SMB, surface mass balance at sea level decreases from an initial value of −0.5 m a−1 m−2 to a final value ranging from −0.9 (SMB_010) up to −4.1 m a−1 m−2 (SMB_090) by the end of the model run, representative of the most negative values of surface mass balance at the margin of the contemporary Greenland Ice Sheet [Ettema et al., 2009]. In the third set of experiments, end‐of‐run rates of submarine melt from 8 (SM_2) to 1000 m a−1 (SM_250) are investigated, providing a range representative of subaqueous melting observed for freshwater‐terminating glaciers [Hochstein et al., 1998] up to the largest order of magnitude in submarine melt rates documented for the present‐day Greenland Ice Sheet [Rignot et al., 2010; Enderlin and Howat, 2013]. Linear and step changes in submarine melting are investigated (Table 2), as described for experiment set 2. The fourth set of experiments combines multiple forcings of crevasse water depth, surface mass balance, and submarine melting. Paired forcing combinations for selected values of linear change are used to investigate the effect of multiple forcings on the rapidity of retreat and to better illustrate the relative sensitivity of ice stream dynamics to each of the paired forcings. Finally, since the above described experiments only address the impact of increased melting on surface and englacial processes, a brief investigation of meltwater‐enhanced basal sliding is also conducted. A triple combination of forcings is applied for this purpose, using low‐to‐moderate end‐of‐run maximum values of crevasse water depth (135 m), surface mass balance (−1.7 m a−1 m−2), and submarine melt (20 m a−1). A basal sliding enhancement was imposed 1000 years into the model runs, reducing the basal friction parameter, μ, in equation 1, by 20% and 30% (CWD10_SMB030_SM5_BF08 and CWD10_SMB030_SM5_BF07, respectively). A control run with full basal friction (CWD10_SMB030_SM5_BF1) is conducted for comparison.

4 Results 4.1 Sensitivity to System Geometry: Sea Level, Topography, and Catchment Width Forcing the model with RSL change alone produced only a very small retreat of the grounding line of approximately 21 km over 2000 years (Figure 2), suggesting that it is unlikely that sea level alone could force grounding line retreat from a stable Younger Dryas position. Accordingly, further sensitivity testing to catchment geometry used surface mass balance run SMB_040 as a control against which to compare modeled response. An initial run in which we remove the RSL forcing further confirms that modeled Bothnian Sea deglaciation is only weakly sensitive to sea level, with the total grounding line retreat distance deviating by only 9.3% between runs with and without RSL change. Grounding line retreat for experiments in which ice stream width was increased and decreased by 25% and 50% deviated very little from control run SMB_040 (Figure 2). The final grounding line positions for these four experiments, at 9.65 ka B.P., ranged between 500 m and 18 km from that of SMB_040, between 2.95 and 0.08% of the total retreat distance. Experiments for decreased width produced the largest, although still small, deviation from the control. These results convey a relative insensitivity to changes in ice stream width within the minimum and maximum ranges as drawn in Figure 1, such that any error in defining the average width for use in subsequent sensitivity tests is assumed to not significantly affect modeled retreat rates. Enhancing tilt of the bed via a 10 and 20% increase in slope in the direction of the ice divide resulted in final grounding line positions exceeding run SMB_040 by 2 km and 57 km, respectively (Figure 2), with considerably less spatial variability preceding the very end of the model run. This amounts to between 0.34 and 9.35% of the total retreat distance, with the grounding line reaching the final position of the surface mass balance‐driven control run approximately 50 and approximately 100 years earlier. In this setting, therefore, the modeled ice stream exhibits relatively little sensitivity to width and bed tilt. We can thus be reasonably confident that the results of our following analyses, which all adopt the control run geometry, will represent the major controls on deglaciation. 4.2 Sensitivity to Atmospheric and Submarine Melting 4.2.1 Individual Forcings The modeled responses of the grounding line position to linear and stepwise changes in crevasse water depth, surface mass balance, and submarine melting are illustrated in Figure 3. The application of sufficiently large individual linear forcings of both crevasse water depth (175 m or more by the end of the run) and surface mass balance (−3.1 m a−1 m−2 or more by the end of the run) result in full retreat of the entire flow line by the end of the model period (Figures 3a and 3c), and all three forcing mechanisms are individually capable of driving retreat of ice out of the present‐day Gulf of Bothnia. While we note that the upper portion of the flow line is not, strictly, appropriate for the known local retreat geometry and timing, we report here on the relative patterns of retreat produced by the different forcing mechanisms over the present‐day marine portion of the flow line, whose geometry is well constrained. Figure 3 Open in figure viewer PowerPoint Grounding line retreat sensitivity to crevasse water depth (CWD), surface mass balance (SMB), and submarine melting (SM) for (a, c, and e) linear and (b, d, and f) stepwise forcing experiments. The black arrows indicate the timing of step changes in forcing. Over the main portion of the trunk of the Bothnian Sea (i.e., approximately 430–700 km along flow line), the modeled retreat rate is rather insensitive to the magnitude of crevasse propagation, when this is forced by a linearly increasing crevasse water depth. The slope of the plotted retreat trajectories in Figure 3a is consistent between model runs. The magnitude of crevassing instead determines how long the grounding line pauses at particular locations and therefore controls the absolute timing of retreat rather than the rate. Higher rates of crevasse propagation succeed in driving retreat of the grounding line to the Härnösand Deep (approximately 400 km along flow line) within our modeled time frame, at which point runaway deglaciation faster than the deglacial chronology of Stroeven et al. [2016] is triggered. Stepped forcing in crevasse propagation, in contrast, produces rather different retreat trajectories depending on the magnitude of the forcing (Figure 3b). In contrast to crevasse propagation, increasing the magnitude of the surface mass balance forcing steadily accelerates the rate of grounding line retreat (Figure 3c). Runs with increasingly negative surface mass balance exhibit less local variability in retreat rate, with fewer and less pronounced inflections in the plotted retreat trajectories under more negative surface mass balance scenarios. Surface mass balance in excess of −1.9 m a−1 m−2 (SMB_035) at the end of run is required to drive the grounding line out of the present‐day Bothnian Sea, while compared to the known deglacial chronology (Figure 1), surface mass balance scenarios of −3.1 to −4.1 m a−1 m−2 (SMB_065 and SMB_090) by the end of the run achieve full deglaciation within an appropriate time frame. Submarine melt‐forced retreat of the ice stream results in relatively little grounding line response to end of run values of up to 100 m a−1 (SM_25; Figure 3e). Only exceptionally high submarine melt rates (200 m a−1 or more) can force complete retreat of this catchment to the terrestrial environment. Figure 4 illustrates the evolution of the ice surface profile and ice flow velocity for selected individual linear forcings of crevasse water depth, surface mass balance, and submarine melting. These three experiments retreat to or close to the prescribed ice divide over the 2000 year model run, with experiment CWD_20, with a maximum water depth of 175 m, fully deglaciating within the final 50 year time step under very rapid retreat. Terminus retreat forced by an increase in crevasse water depth is characterized by two periods of relative retreat slowdown (Figure 4a) coinciding with a reduction in ice surface velocities (Figure 4d), both of which fall on a gentle reverse bed slope just prior to its abrupt steepening. The first occurs at approximately 720 km from the ice divide and the second at approximately 450 km from the ice divide. This pronounced deceleration in the rate of grounding line retreat coincides with a narrowing of the ice stream width, but is seen only in runs for which crevasse water depth is the sole, or one of two applied forcings (Figures 3a, 3b, 5a, and 5c), and not when the model is forced by a change in submarine melting or surface mass balance alone. Immediately preceding full deglaciation in the final time step of run CWD_20 there is a large peak in ice surface velocity, reaching approximately 6000 m a−1 (Figure 4d). This rapid retreat in the final 50 years could be attributable to a debuttressing effect caused by loss of the ice shelf in the final time step, to a reduction in basal and/or lateral drag, or to acceleration of ice surface velocities in response to a steepening ice surface profile. When forced by decreasing surface mass balance alone, the ice stream surface profile lowers at a considerably greater rate than when forced by crevasse water depth or submarine melting (Figure 4b). The decreasing mass balance and thus ice flux leads to a reduction in ice surface velocities along the flow line, and without additional increases in calving or submarine melt, a large floating ice tongue can form, stretching up to approximately 100 km in length. Grounding line retreat forced by submarine melt alone (Figure 4c) is steady in comparison to crevasse water depth‐forced retreat, with less variation in ice surface velocities (Figure 4f), and no major periods of retreat slowdown. Figure 4 Open in figure viewer PowerPoint Comparison of terminus retreat and ice surface evolution for experiments: (a) CWD_20, (b) SMB_040, and (c) SM_50, plotted for each 50 year time step over a total of 40 time steps (2000 years). The thickness of the bed (black line) depicts the total RSL change over the model run. (d–f) The ice surface velocities along the flow line at each 50 year time step for experiments CWD_20, SMB_040, and SM_50, respectively. Note that the experiment CWD_20 retreats fully to the ice divide within the final 50 year time step. Figure 5 Open in figure viewer PowerPoint Grounding line retreat for combined linear forcing experiments, testing sensitivity to (a) crevasse water depth (CWD) and submarine melting (SM), (b) surface mass balance (SMB) and submarine melting (SM), and (c) crevasse water depth (CWD) and surface mass balance (SMB). The application of linearly increasing forcings results in retreat rates which are nonlinear. This indicates sensitivity of grounding line retreat to local changes in basal topography. There are consistent inflections in the grounding line responses of many of the model runs, including at approximately 780, 745, 720, 625, and 485 km along the flow line, which correspond to abrupt changes in local basal slope (Figure 3). These bed topography changes therefore appear to control small‐scale grounding line responses regardless of the mode of “external” applied forcing. Stepwise forcing experiments, in which larger increases were imposed every 500 years to investigate response to abrupt forcings, yield a grounding line response which only occasionally responds in a stepwise fashion. Response to step changes is most pronounced for crevasse water depth experiments CWD_step20 and CWD_step30 (Figure 3b), where acceleration in grounding line retreat is enhanced following step increases at 10.6 and 10.1 ka B.P. Local variation in retreat rates is subdued during these step‐driven periods of rapid retreat, suggesting that while the modeled ice stream is adapting to these sudden changes, the influence of local topography is overridden by other controls. Neither surface mass balance nor submarine melting rates drive any abrupt retreat events at the times when an abrupt (stepped) forcing is applied. In the case of the former this may be due to surface mass balance having a time‐integrated effect; thus, an abrupt reduction in surface mass balance may not have an instantaneous effect on the grounding line position. The above results indicate that individual forcings, if they are sufficiently large, can lead to full retreat by the end of the model period. These results also demonstrate that while sufficient step changes in crevasse propagation via increasing crevasse water storage forces acceleration of retreat rates, this response is short‐lived. Within our suite of model experiments, step changes in individually applied forcings do not, alone, cause a catastrophic deglaciation event. 4.2.2 Paired Forcings Each forcing mechanism can drive complete retreat of the modeled Bothnian Sea flow path if its magnitude is sufficiently high. It is instructive, therefore, to examine paired forcing combinations to allow the relative influence of crevasse water depth, surface mass balance, and submarine melting on retreat trajectories to be evaluated more closely. Figure 5 illustrates the grounding line retreat for each combination of forcings, where forcings are kept low to moderate in magnitude (Table 2) in order to allow us to best judge sensitivity without one signal artificially overriding another. Individually, none of these modest forcings would cause retreat to more than approximately 500 km from the ice divide after 2000 years (central Bothnian Sea), but in certain combinations they amount to full or close‐to‐full retreat. The most striking result of pairing forcings is that the effects of submarine melt are negligible. Figure 5a shows that crevasse water depth exerts a substantially greater control on modeled grounding line retreat than submarine melting. Grounding line retreat trajectories cluster tightly according to the magnitude of crevasse water depth, while varying the submarine melt rate has very little effect. Only experiment CWD15_SM10 exhibits notable retreat of the final grounding line position in response to increased submarine melting, and this is clearly coincident with a steep reverse bed slope. While runs pairing surface mass balance and submarine melt do not exhibit the same tight clustering of retreat trajectories, these responses plot in an order dictated by surface mass balance magnitudes (Figure 5b), and therefore, we interpret the surface mass balance forcing to similarly dominate over submarine melting. In contrast, it is difficult to separate the relative influence of crevasse water depth and surface mass balance (Figure 5c). Grounding line retreat trajectories crosscut one another, and neither mechanism appears to dominate. We note that our greatest (most negative) surface mass balance forcing of −0.3 m a−1 m−2 over 500 years (−1.7 m a−1 m−2 by end of run) drives more rapid and extensive retreat than an increase in crevasse water depth in combination with smaller changes in surface mass balance. It is apparent, however, that within the magnitudes of forcings explored here, surface melt‐related processes dominate over submarine melting in driving grounding line retreat in the Bothnian Sea. 4.2.3 Full Forcing and Sensitivity to Basal Friction While enhanced surface melting appears to be an important driver of ice retreat in the Bothnian Basin, we have thus far only considered its effect via surface mass balance reduction and surface crevasse propagation. Surface melting is not physically coupled within the model to any basal‐sliding response mechanism. Reducing basal friction allows for a simple investigation of the possible influence of increased surface meltwater production and delivery to the subglacial domain. Three final multiple forcing experiments were conducted with a fully combined (triple) forcing setup, comprising a control run with baseline basal friction, and two runs with reduced basal friction (Figure 6). With the baseline basal friction (i.e., as in all previous simulations), the combination of forcings yields a pattern of margin retreat that sticks persistently over the Åland sill (approximately 720 km), pauses again at approximately 600–620 km, and thereafter the retreat rate increases. The accompanying ice velocities, which peak initially approximately 670 km from the ice divide immediately following rapid retreat over the deep water of the Åland Sea, are reduced to moderate‐fast flow speeds through the central Bothnian Sea and rise again as the retreat rate increases over the reverse bed and deep basin of Härnösandsdjupet (approximately 500–350 km; Figures 6a and 6d). We find that a 20 and 30% reduction in basal friction does not significantly alter either the spatial pattern of retreat or the magnitude of grounding line retreat accomplished in 2000 years, though increases some minor spatial variability (i.e., standstills and retreat steps). Absolute ice flow velocities increase with the reduction in basal friction, particularly between approximately 600 and 670 km, but the spatial pattern of velocity evolution is similarly insensitive to basal friction. Figure 6 Open in figure viewer PowerPoint Comparison of terminus retreat and ice surface evolution for enhanced basal sliding experiments: (a) CWD10_SMB030_SM5_BF1, (b) CWD10_SMB030_SM5_BF08, and (c) CWD10_SMB030_SM5_BF07, plotted for each 50 year time step over a total of 40 time steps (2000 years). (d–f) The ice surface velocities along the flow line at each 50 year time step for experiments CWD10_SMB030_SM5_BF1, CWD10_SMB030_SM5_BF08, and CWD10_SMB030_SM5_BF07, respectively. (g) The locations of ice stream termination, sticky zones, and basal crevassing, as inferred from geomorphological analysis of multibeam data.

5 Discussion The results described above indicate that each individual mechanism could lead to complete or close to complete retreat if the applied forcing is sufficient. However, it is likely that such values may become unreasonable in order to achieve this. We designed the range of experimental values according to those reported in both/either contemporary or paleo settings, specifically to meet our sensitivity testing objectives, but it does not follow that this range should be applicable to the case of the Bothnian Sea. This applies particularly to submarine melting, for which even a melt rate of 1000 m a−1, albeit through a very simple implementation, does not result in full retreat of the ice stream (Figure 3e). While such magnitude of submarine melting may be applicable to particular individual outlets of the present‐day Greenland Ice Sheet [Rignot et al., 2010; Enderlin and Howat, 2013], it is an exceptionally high rate compared to other Greenland and Antarctic catchments [Jacobs et al., 2011; Enderlin and Howat, 2013; Rignot et al., 2013], and it is highly unlikely that melting of this extent was active in the Gulf of Bothnia given that climate was likely still cooler than present during deglaciation. During the time period of Bothnian Sea deglaciation (approximately 11.6–10.3 ka), exchange of Baltic‐Bothnian waters with the North Sea was extremely limited, confined to the ever‐shallowing and ever‐narrowing Närke Strait across central Sweden as land uplift proceeded to close the basin outlet (Yoldia Sea to Ancylus Lake transition [Björk, 1995; Andrén et al., 2011]). Limited marine inflow created a brackish environment for a maximum of approximately 300 years, 11.4–11.1 ka [Andrén et al., 2011]; otherwise, the Bothnian Sea was a freshwater body fed by cold glacial meltwater. Without proglacial circulation driven by either water mass exchange or plume‐driven convection [Jenkins, 2011], we expect low transfer of heat from the proglacial water body to the ice front and therefore very low melt rates. In the range of values that are realistic in this setting, we conclude from our model responses that it is likely that the Bothnian Sea Ice Stream was relatively insensitive to submarine melting. Through quantitative analysis of mass loss partitioning for the triple forcing experiment (Figure 6a), it is clear that the contribution of submarine melting to overall mass loss is relatively small in comparison to surface mass balance and calving (Figure 7). However, while absolute contribution to mass loss is low, small peaks in submarine melting precede some large calving‐driven retreat episodes in this experiment, suggesting that it may have a more important role in triggering retreat from topographic pinning points. Figure 7 Open in figure viewer PowerPoint Rate of mass loss at each 50 year model time step for experiment CWD10_SMB030_SM5_BF1, partitioned by calving, surface mass balance, and submarine melt. Note that calving is a by‐product of crevasse water depth but that calving events are not solely attributable to crevasse water depth. Our sensitivity tests to bed geometry revealed that retreat of this ice sheet sector was likely relatively insensitive to RSL change and associated isostatic adjustment. This insensitivity of modeled grounding line retreat to either submarine melting or RSL change indicates that this marine ice sheet sector was governed only minimally by direct marine processes. Conversely, the modeled grounding line retreat of the Bothnian Sea Ice Stream is shown to be highly sensitive to surface melting, which acts both through enhanced crevasse propagation due to increased water infilling of surface crevasses (i.e., hydrofracture) and through an increasingly negative surface mass balance and associated thinning of the ice surface profile. We cannot presently separate the influence of surface mass balance and crevasse propagation. Qualitatively, crevasse water depth appears to affect the timing of retreat more than the rate via control on the duration of episodes, where the grounding line appears to undergo limited movement. Surface mass balance affects retreat rates in a steadier fashion. A large negative surface mass balance (end‐of‐run values of −3.1 to −4.1 m a−1 m−2) is independently capable of driving complete Bothnian Sea retreat within a time frame matching the terrestrial clay varve chronology for deglaciation [Stroeven et al., 2016]. Such surface mass balance rates are in fact rather moderate compared to an estimated −5 to −9 m a−1 m−2 for the southern Laurentide Ice Sheet during early Holocene [Carlson et al., 2009], and the geomorphological record exhibits abundant evidence of high‐discharge, well‐connected subglacial meltwater drainage [Greenwood et al., 2016]. Implementation of stepwise reductions in surface mass balance, however, did not drive instantaneous response or catastrophic retreat, possibly due to the time‐integrated effect of surface mass balance (Figure 3d). Widespread basal crevasse squeeze ridges additionally point to a highly fractured ice body, likely vulnerable to surface‐melt enhanced hydrofracture. While we cannot yet quantitatively disentangle mass balance and hydrofracture effects on retreat rate, our model implicates ice surface melting in driving deglaciation, as illustrated by its contribution to overall mass loss (Figure 7). Increased melting and increased crevasse propagation are closely linked on contemporary ice sheets [Clason et al., 2015], and collectively act as triggers for secondary forcings including increased calving and ice cliff failure [Benn et al., 2007; Pollard et al., 2015; Colgan et al., 2016], and changes in efficiency of the subglacial hydrological system [Cowton et al., 2013; Meyaud et al., 2014]. The extent to which increased surface‐to‐bed meltwater transfer influences annual ice surface velocities remains contentious, and the extent to which increased efficiency of the subglacial drainage system may offset enhanced basal lubrication on the Greenland Ice Sheet remains uncertain [Sundal et al., 2011; Sole et al., 2013; van de Wal et al., 2015]. Given the likely importance of meltwater in the Bothnian Sea system, we employed a simple proxy for meltwater‐enhanced sliding to allow for a first‐order exploration of the possible effects on ice dynamics and retreat. While reducing basal friction did result in increased velocities along the flow line overall, it did not produce any significant variation in the spatial pattern of velocities along the flow line. Additionally, we do not observe any differences in the overall duration of deglaciation, despite higher velocities accompanying low friction; the velocity response, in itself, does not make the system any more vulnerable to rapid retreat. In this setting, enhanced flow rates are offset by an enhanced duration of grounding line pinning at those localities prone to grounding line retreat slowdown (Figure 6). Consistent areas of grounding line retreat slowdown and of increased retreat rates are produced across the suite of model experiments. We interpret this pattern to indicate model sensitivity to local variations in basal topography and ice stream width, despite the relative insensitivity to the absolute magnitude of these factors at the catchment scale (Figure 2). Under almost all forcings applied in this study the margin retreat slows approximately 720 km along the flow line for a short period of approximately 250 years, immediately preceding rapid retreat of approximately 50 km within 50 years through an area of deep water in the Sea of Åland. Grounding line retreat slowdown also occurs at approximately 450 km along flow, on the gentle reverse slope distal to the Härnösand Deep, most commonly and markedly when crevasse propagation is part of the applied forcing. This is clearly the case in both single and paired forcing runs (Figures 3-5), although the full forcing scenario only displays a brief (approximately 50 year) slowdown of retreat under the enhanced basal sliding scenario in Figure 6c. In the full forcing run and in the paired CWD‐SMB forcing, there is a phase of markedly reduced retreat rate approximately 650–600 km along flow line. The geomorphological and geological records of the Bothnian Sea identify (i) a paleo‐ice stream pathway, with high ice flow velocities limited to a stretch approximately 250–620 km along our modeled flow line; (ii) the terminal zone of the active ice stream, located in the central‐southern Bothnian Sea approximately 600–620 km along flow line (Figure 6g); (iii) a zone of relatively enhanced basal stability immediately south of the Härnösand Deep (Figure 1); and (iv) a 5–30 m thick sequence of glacial sediments in the main trunk of the Bothnian Sea, while the distal end of our flow line (from approximately 660 km) passes over crystalline bedrock with only a thin (<10 m) sediment cover. The modeled positions of slowdowns in grounding line retreat, and the general lack of evidence for prolonged periods of retreat slowdown through the Bothnian Sea, are consistent with these properties. The first modeled position of minor slowdown, across the Åland shallows, occurs on the crystalline substrate where increased bed roughness may be expected to locally reduce ice flow velocities and stabilize the grounding line [Rippin et al., 2011; Livingstone et al., 2012] and where small glacial lineations are only weakly developed. Additionally, this position coincides with a narrowing of the ice stream width (Figure 1), which is reported elsewhere to permit a slowdown of retreat even on reverse bed slopes [Jamieson et al., 2012; Gudmundsson et al., 2012]. The island of Åland itself has, in the past, been proposed as a “sticky spot” [Holmlund and Fastook, 1993], acting to inhibit ice flux through the south‐central sector of the FIS, which has led to overly rapid retreat in previous modeling studies [Clason et al., 2014]. Slowdown of retreat approximately 450 km along flow line, just south of Härnösand Deep, emerges consistently from our model runs and is also indicated by the appearance of transverse wedge‐like ridges superimposed on MSGLs [Greenwood et al., 2016] and which record a zone of reduced flow and margin stability after a period of ice streaming. Since modeled grounding line retreat only slows across this zone when surface crevasse propagation is a component of the forcing, we implicate enhanced crevassing as a control on stability in this system. Deeper crevasses lead to an increase in calving, maintaining a steep ice sheet profile and restricting formation of a floating tongue, resulting in a grounding line that is more resistant to buoyancy and thus potentially more likely to stabilize over reverse bed slopes. With paired atmospheric forcings (surface mass balance and crevasse water depth) and with a fully combined (triple) forcing, but notably in none of the individual forcings nor the alternative pairs, the model predicts a prominent zone of grounding line retreat slowdown at approximately 650–600 km along flow line, which lasts for 200–250 years. This is consistent with the geomorphological record, which reveals the terminal zone of the Bothnian Sea Ice Stream at approximately 620–600 km (Figure 6g). In this zone, sets of MSGLs crosscut subparallel to one another, in contrast to further up the ice stream trunk where lineations display a high degree of parallel conformity (e.g., Figure 1d). While there is no evidence for a grounding zone wedge or other ice‐marginal deposits here, the lineation overprinting requires that there was a persistent ice stream terminus at approximately 620–600 km, maintaining its position for sufficient time to allow the marginal flow direction to shift back and forth a few degrees. Our modeling supports this central‐southern Bothnian Sea position of a moderately stable grounding line, and its appearance in only the atmospherically forced runs again attests to the likely importance of ice surface processes in governing this ice sheet sector. The patterns of ice velocity evolution, on the other hand, are an imperfect match to the reconstruction of ice streaming from the glacial landform record. That there is only an MSGL record of ice streaming up to 620 km along our flow line, beyond which lies a separate group of small (often bedrock) drumlins and crag and tails with a flow direction offset from the MSGLs, suggests that ice streaming occurred as a distinct event once retreat into the Bothnian Sea had already begun. Furthermore, this ice stream event did not extend the full length of the Bothnian Sea [Greenwood et al., 2016]. We would therefore expect to see initially low flow velocities, with an abrupt shift to high velocities once the margin reached approximately 620 km. Modeled flow velocities, in contrast, peak immediately upon withdrawal from the Åland Sea, at approximately 680–690 km, and then decline toward 600 km, although still at rates typical of ice streaming, i.e., 500–1000 m a−1. Modeled velocities peak again around 500 km along flow line, in the central Bothnian Sea, which is coincident with where the length of MSGLs reach their peak. In its present form the model does not capture the Bothnian Sea Ice Stream as an “event” which abruptly activates during the Bothnian Sea deglaciation or the precise timing of southern FIS deglaciation known from the terrestrial clay varve chronologies (summarized in Stroeven et al. [2016]). However, our present implementation of atmospheric and marine drivers is designed only to test sensitivity to retreat triggers, and we conclude that it is likely that the ice streaming event did not arise simply as a nonlinear response to steady climate amelioration. Our model experiments highlight surface melting as an important control on Bothnian Sea grounding line retreat. Despite limitations associated with a one‐dimensional model, our results shed new light upon the mode of marine‐based ice sheet retreat in the southern FIS and, by extension, other low‐relief terrains associated with large, proglacial water bodies. That the Bothnian Sea Ice Stream is a calving ice stream that most likely collapsed without significant sensitivity to marine forcing (i.e., sea level change and submarine melting) highlights the need to considerably improve our understanding of how subaqueous ice sheet sectors respond to both atmospheric and marine forcings.

6 Conclusions We applied a well‐established flow line model to the case of the Bothnian Sea Ice Stream, validated against geomorphological evidence that has been mapped and analyzed from recently collected high‐resolution bathymetric data. Results of sensitivity testing to both single and multiple forcings of crevasse water depth, surface mass balance, and submarine melting point toward retreat governed largely by increased surface meltwater production in response to climatic changes following the Younger Dryas. Submarine melting and relative sea level change are found to have only minimal effect on grounding line retreat in our simulation. We thus highlight the Bothnian Sea as a case where despite the marine setting, retreat of this large ice sheet sector was likely governed primarily by atmospheric, rather than marine processes. While the modeling presented here applies very simple implementations of calving, surface, and submarine melting, our results stress the importance of better understanding the multitude of responses to climatic and environmental controls exhibited by marine‐terminating ice sheets. Our suite of model experiments indicates no major standstills after the Younger Dryas but short‐lived periods of grounding line retreat slowdown and of accelerated retreat. The locations of slowdowns are consistent with the geomorphological record and are interpreted to indicate sensitivity to local variations in basal topography and in ice stream width. While the relative influence of meltwater‐forced crevasse propagation and surface mass balance cannot yet be separated, we find that the modeled fit with the geomorphological record of grounding line stabilization is strongest when both processes are fully incorporated. The model does not replicate the activation of the Bothnian Sea Ice Stream as a discrete, fast‐flow event that is triggered once retreat is already underway. Since activation does not arise from either linearly or regularly stepped increases in forcing, we infer that the ice stream may have switched on in response to a specific event forcing not accounted for in our current modeling setup. Implementation of a realistic paleo‐climate forcing is necessary to evaluate the temporal and spatial responses of the Bothnian Sea Ice Stream to climatic variations under a period of rapid environmental change. The availability of high‐resolution topographic and bathymetric data in this region, in concert with geomorphologically validated numerical modeling, provides a unique opportunity to improve our understanding of ice stream dynamics and collapse in a major basin of the Fennoscandian Ice Sheet, until now a gap in the deglacial history of the region.

Acknowledgments We gratefully acknowledge the financial support of the Swedish Radiation Safety Authority (P.H.), the Carl Trygger Foundation for Scientific Research (C.C.C.), the Geological Survey of Sweden, and the Swedish Research Council (S.L.G.). The contribution of J.M.L. was partly supported by FORMAS grant 2013‐1600. The Swedish Maritime Administration and the Geological Survey of Sweden kindly provided access to geomorphological and geological data. Bathymetric and topographic data can be accessed via the Baltic Sea Bathymetry Database (http://www.bshc.pro/) and Lantmäteriet (http://www.lantmateriet.se/), respectively. We wish to thank Francesco Muschitiello and Arjen Stroeven for their discussions which helped to improve the study. We also thank Nick Golledge, Olga Sergienko, and an anonymous reviewer for their helpful comments.