The influence of glacial landscape evolution on Scandinavian ice-sheet dynamics and dimensions

The Scandinavian topography and bathymetry have been shaped by ice through numerous glacial cycles in the Quaternary. In this study, we investigate how the changing morphology has influenced the Scandinavian ice sheet (SIS) in return. We use a higher-order ice-sheet model to simulate the SIS through a glacial period on three different topographies, representing different stages of glacial landscape evolution in the Quaternary. By forcing the three experiments with the same climate conditions, we isolate the effects of a changing landscape morphology on the evolution and dynamics of the ice sheet. We find that early Quaternary glaciations in Scandinavia were limited in extent and volume by the pre-glacial bathymetry until glacial deposits filled depressions in the North Sea and built out the Norwegian shelf. From middle–late Quaternary (∼0.5 Ma) the bathymetry was sufficiently filled to allow for a faster southward expansion of the ice sheet causing a relative increase in ice-sheet volume and extent. Furthermore, we show that the formation of The Norwegian Channel during recent glacial periods restricted southward ice-sheet expansion, only allowing for the ice sheet to advance into the southern North Sea close to glacial maxima. Finally, our experiments indicate that different stretches of The Norwegian Channel may have formed in distinct stages during glacial periods since ∼0.5 Ma. These results highlight the importance of accounting for changes in landscape morphology through time when inferring ice-sheet history from ice-volume proxies and when interpreting climate variability from past ice-sheet extents.

1 Introduction

Ice holds the power to transform landscapes and constituted a major geomorphological agent in northern Europe during the Quaternary (last 2.6 Ma) where recurring glacial cycles shaped the present-day landscape. Indeed, the topography and bathymetry in and around northern Europe reveal the extensive impact of its rich glacial history, with deep fjords and U-shaped valleys attesting to the accumulated effect of widespread glacial erosion and terminal moraines indicating the extent of past ice sheets (Hughes et al., 2016; Stroeven at al., 2016). The Eurasian ice-sheet complex covered much of the British Isles, all of Scandinavia, and much of northern Europe including parts of Germany, Poland, Russia, and the Baltic through multiple glacial cycles since 1 Ma (Batchelor et al., 2019). During the Last Glacial Maximum (LGM), the ice-sheet complex, consisting of the Scandinavian ice sheet (SIS), the Barents Sea ice sheet (BSIS), and the British–Irish ice sheet (BIIS), contained an ice volume corresponding to $\sim \mathrm{18.4}\pm \mathrm{4.9}$ m sea-level equivalent (Simms et al., 2019). On a global scale, the pace of these glacial cycles results from solar insulation variations combined with feedback mechanisms and internal dynamic effects in the climate system, in part caused by the ice sheets themselves (Hughes and Gibbard, 2018). Differences in both ice volume and the extent of ice sheets between glacial cycles (Fig. 1) can also be attributed to variations in moisture supply through complex global atmosphere–ocean–ice interactions (e.g., Batchelor et al., 2019; Hughes and Gibbard, 2018), with topography and proximity to the ocean being key factors determining the spatial distribution of moisture to an ice sheet. Studies on glacial landscape evolution have indicated that glacial erosion and deposition can also influence ice-sheet dynamics as well as ice volumes and extent (e.g., Kessler et al., 2008; Kaplan et al., 2009; MacGregor et al., 2009; Egholm et al., 2009, 2012a, b, 2017; Anderson et al., 2012; Pedersen and Egholm, 2013; Pedersen et al., 2014; Clague et al., 2020; Mas e Braga et al., 2023). But until now, these studies have been limited to synthetic landscapes and/or limited spatial scales (smaller glaciers and ice caps). A few ice-sheet scale models are starting to consider glacial erosion (e.g., Patton et al., 2022), but the effects of long-term Quaternary landscape evolution on ice-sheet dynamics are still to be explored on a large scale for realistic landscapes and ice-sheet configurations. Understanding the influence of landscape evolution on ice-sheet dynamics requires the reconstruction of landscapes that existed prior to or at earlier stages of glacial erosion, something that can be approached using source-to-sink studies, utilizing offshore sediment volumes of a glacial origin (e.g., Steer et al., 2012; Paxman et al., 2019; Pedersen et al., 2021).

Figure 1 Overview map of the model domain. The maximum plausible extent of the Fennoscandian ice-sheet complex during the Last Glacial Maximum (LGM; black line) and the Penultimate Glacial Maximum (MIS 6; dashed red line) are overlaid (Batchelor et al., 2019). Also shown is the approximate location of the LGM ice-divide position (Olsen et al., 2013).

In this work, we focus on the well-studied Scandinavian region and investigate how the SIS may have changed its behavior because of Quaternary landscape evolution. We use a higher-order ice-sheet model to investigate how large-scale glacial morphological features have influenced the development and dynamics of the SIS over a glacial cycle at two key times during the Quaternary: (1) before the inception of major glaciations in the beginning of the Quaternary (PREQ ∼2.6 Ma) and (2) during the middle–late Quaternary (MLQ ∼0.5 Ma) where major pre-glacial features in the bathymetry around Scandinavia had been filled with glacial deposits (Dowdeswell and Ottesen, 2013). Importantly, we do not intend to reconstruct realistic SIS configurations for these past time periods, but rather keep the climate forcing consistent between experiments, in order to isolate how changes in bed morphology has impacted SIS dynamics and extent. This allows us to (i) explore how morphological changes can influence the dynamics, extent, and volume of the ice sheet, independent of the climatic forcing, and (ii) gain insight into how ice-volume proxies could be influenced by glacial landscape evolution.

For the early Quaternary, we adopt the pre-glacial landscape reconstructions provided for the Scandinavian region by Pedersen et al. (2021) that include (i) the absence of glacially generated sediments offshore, (ii) infill of over-deepened fjords and glacial valleys onshore, (iii) a reconstructed wedge of older Mesozoic and Cenozoic sediments on the inner shelf that is assumed to have been eroded by glacial activity within the Quaternary (e.g., Hall et al., 2013), and finally, (iv) adjustments of the landscape owing to erosion- and deposition-driven isostatic changes and dynamic topography (Pedersen et al., 2016).

In addition to this pre-glacial reconstruction, which explores an entirely different offshore bathymetry and onshore Scandinavian landscape, we also consider the more subtle effects of the large glacial troughs that have been carved into the shelf bathymetry by ice streams since the middle–late Quaternary. One of the most notable of these glacio-morphological features offshore Scandinavia is The Norwegian Channel (Fig. 1). This channel is believed to have been formed by ice-stream activity sometime since 1.1 Ma (e.g., Sejrup et al., 2003), with studies suggesting that ∼90 % of the deposits funneled through the channel and into the North Sea Fan were deposited within the last ∼0.5 Ma (Hjelstuen et al., 2012). Recently, it has been argued that the channel formed before ∼0.35 Ma (Løseth et al., 2022). An erosional unconformity at the base of the channel is draped by post-LGM sediments, suggesting that the channel experienced erosion within the last glacial cycle (Hjelstuen et al., 2012).

2 Methods

For the numerical experiments presented in this study, we use the depth-integrated second-order shallow-ice approximation iSOSIA (Egholm et al., 2011, 2012a, b). We conduct our experiments by simulating a full glacial cycle of 120 ka on different topographies. In the following section we will present the numerical model, the model setup, and the experimental design.

2.1 Modeling the Scandinavian ice sheet

The ice flow in iSOSIA is governed by a second-order approximation of the equations for Stokes flow (e.g., Egholm et al., 2011). The velocities are depth integrated to yield a 2D one-layer ice model, implemented here using a regular grid (e.g., Egholm and Nielsen, 2010). The second-order nature of the approximation ensures that ice velocities depend non-linearly on ice thickness, ice-surface gradients, as well as longitudinal and transversal horizontal stress gradients (Egholm et al., 2011, 2012b). Details on the iSOSIA model, including the importance of the higher-order ice dynamics involved, have been described in depth elsewhere (Egholm and Nielsen, 2010; Egholm et al., 2011, 2012a, b).

The depth-integrated ice-creep velocity is calculated using temperature-dependent Glen's flow with a stress exponent, n, equal to 3:

$$\begin{array}{}\text{(1)}& \dot{{\mathit{\epsilon }}_{ij}}=A_{\mathrm{flow}}{\mathit{\tau }}_{\mathrm{e}}^{n-\mathrm{1}}{s}_{ij},\end{array}$$

where $\dot{\mathit{\epsilon }}$ is the strain rate tensor, ij denotes the components of the tensor, A_flow is the ice flow parameter, τ_e is the effective stress, and s is the deviatoric stress tensor (Egholm et al., 2011). The ice flow parameter A_flow is dependent on the depth-averaged temperature of the ice using an exponential relationship:

$$\begin{array}{}\text{(2)}& A_{\mathrm{flow}}=A_{\mathrm{0}}\exp \left({\displaystyle \frac{-Q}{RT}}\right),\end{array}$$

where A₀ is a flow constant, Q is an activation energy, R is the gas constant, and T is the temperature relative to the pressure melting point (e.g., Zeitz et al., 2020). A₀ and Q have different values above and below $T=-\mathrm{10}$ °C (see Table 1). A simple Weertman sliding scheme is used to calculate the contribution of basal sliding to depth-integrated ice velocities:

$$\begin{array}{}\text{(3)}& u_{\mathrm{b}}=A_{\mathrm{sliding}}{\displaystyle \frac{t_{\mathrm{s}}^{\mathrm{3}}}{N}},\end{array}$$

where u_b is the basal velocity, A_sliding is an ice sliding coefficient, t_s is the bed-parallel shear stress, and N is the effective pressure at the base (Egholm et al., 2011). A_sliding is chosen to give realistic sliding in the order of several hundred meters per year, for example, in fjords or near the shelf edge in the Norwegian Sea, similar to surface velocities in comparable areas of modern-day ice bodies (e.g., Millan et al., 2022). To allow for faster ice flow for soft-bed subglacial conditions (e.g., Gladstone et al., 2020; Han et al., 2021), A_sliding is enhanced by a factor of 5 in offshore regions and onshore in northern Europe where thick, soft sediments cover the bed.

In this study, we focus on grounded ice only, as ice-shelf dynamics are computationally expensive to resolve on the timescales of our experiment and because constraints on ice-shelf extent in middle or early Quaternary glaciations are sparse due to a lack of reliable dates on submarine landforms (e.g., Jakobsson et al., 2016). Some older studies suggest that an ice shelf was present during recent glaciations in the North Atlantic and Arctic regions (Hughes et al., 1977; Lindstrom et al., 1986). However, while ice-shelf stability is sensitive to bathymetric configurations (Bart et al. 2016) and is a deciding factor in grounding line migration, we limit our focus here to large-scale morphological features, such as the Norwegian Channel, created by an ice stream in contact with the seabed (Sejrup et al., 2016). Consequently, we do not consider floating ice in our simulations and remove floating ice by introducing a fast melt rate for ice that does not meet the grounding criterion:

$$\begin{array}{}\text{(4)}& H_{\mathrm{ice}}>\left(\text{SL}+H_{\mathrm{ice}}\right){\displaystyle \frac{{\mathit{\rho }}_{\mathrm{water}}}{{\mathit{\rho }}_{\mathrm{ice}}}},\end{array}$$

where H_ice is ice thickness, SL is local sea level, and ρ_water and ρ_ice are the densities of water and ice, respectively. Mean sea level in the model is varied between the interglacial and glacial maximum (−130 m) using the normalized LR04 Benthic Stack (Lisiecki and Raymo, 2005) as a glacial index. Special boundary conditions are employed at the approximate locations where the SIS meets the BSIS and BIIS by introducing an “ice wall” where the ice flux is zero to emulate divergent ice flow when these ice sheets merge during glacial maxima. At the edges of the model domain, we employ open boundary conditions to allow ice to flow out of the domain. Common model parameters are presented in Table 1.

Table 1 Common parameters in the ice-sheet model and mass balance scheme. Numbers in brackets denote minimum and maximum values.

2.1.1 Mass balance

In the simulations we present here, we assume that the mass balance (${\dot{M}}_{\mathrm{ice}}$) of the ice sheet can be approximated using three components:

$$\begin{array}{}\text{(5)}& {\dot{M}}_{\mathrm{ice}}={\dot{m}}_{\mathrm{acc}}-{\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}},\end{array}$$

where ${\dot{m}}_{\mathrm{acc}}$ is the rate of accumulation, ${\dot{m}}_{\mathrm{s}}$ is the surface melt rate, and ${\dot{m}}_{\mathrm{b}}$ is the basal melt rate (Egholm et al., 2012b). We use a positive-degree-day (PDD) model to estimate accumulation rate and surface melt rate as a function of mean annual temperature, annual temperature variation, and mean annual precipitation at every point in our model domain for every time step (e.g., Magrani et al., 2022).

The yearly temperature variation in each cell is approximated by a sine function based on the mean annual temperature and annual temperature amplitude (see below). The melt rate (in m yr⁻¹) is calculated in the PDD model as

$$\begin{array}{}\text{(6)}& {\dot{m}}_{\mathrm{s}}=m_{\mathrm{PDD}}{T}_{\mathrm{positive}},\end{array}$$

where m_PDD is the positive-degree-day factor multiplied by the sum of positive degrees T_positive each year. Here, we consider a single melting degree factor for both ice and snow since all precipitation is turned into ice after accumulation (based on yearly average rates). The accumulation rate is approximated by

$$\begin{array}{}\text{(7)}& {\dot{m}}_{\mathrm{acc}}={\displaystyle \frac{n_{\mathrm{frost}}}{\mathrm{365}}}\cdot P,\end{array}$$

where n_frost is the number of days with negative temperatures in a year and P is the annual precipitation. The temperature forcing that drives spatial and temporal changes in mass balance in our simulations is based on mean temperature, annual temperature amplitude, and lapse rate, all of which vary across the model domain using spatial gradients that vary in time. Two climate states, a glacial maximum state and an interglacial state, are chosen to represent the extremes of our model, and the spatial gradients of the full glacial cycle of our model simulations are subsequently defined to vary in between these extremes using a glacial index that resembles the normalized LR04 Benthic Stack (Lisiecki and Raymo, 2005) with the glacial maximum in this climate forcing occurring at 18 ka BP. Here, we define spatial (x, y, z) gradients at the glacial maximum using multiple linear regression on MPI-ESM climate model outputs (LGM experiment; Jungclaus et al., 2019). For the interglacial state we define spatial gradients using the ERA-Interim reanalysis data for modern day (Dee and National Center for Atmospheric Research Staff, 2022). Finally, the lapse rate was found to be close to constant, so we keep this fixed at 6.5 °C km⁻¹. With this approach, the temporally varying temperature forcing of the entire grid can be defined from a single grid cell in the lower left corner while still capturing a coastal–continental (east–west) gradient, a polar gradient (south–north), and an altitudinal gradient (lapse rate) in temperature. However, we cannot capture local effects that arise from changes in complex atmospheric circulations patterns over time that might have important implications for glacial dynamics and ice extent (e.g., Liakka et al., 2016; Hughes and Gibbard, 2018).

To represent precipitation in our simulations, we use a climate-corrected modern-day mean precipitation field (Pendergrass et al., 2022), modulating the local precipitation in every grid cell using the following equation:

$$\begin{array}{}\text{(8)}& P=P_{\mathrm{0}}\cdot e^{kTp\cdot \mathrm{\Delta }T},\end{array}$$

where P₀ is the local modern-day (interglacial) precipitation, ΔT is the change in temperature in a cell from the previous time step, and kTp represents the rate of change in precipitation for a change in temperature with a value of 0.029 °C⁻¹. The value of kTp is found by optimization through a comparison between mean precipitation at the LGM in MPI-ESM and mean precipitation in the modern-day ERA-Interim data set. By scaling the precipitation with changes in temperature we can capture some of the effects an ice sheet will impose on moisture supply, by limiting snowfall in the central parts of the ice sheet (Fig. 3d).

Basal melt rate is calculated as the difference between geothermal heat flux from the bed q_b and the heat flux from the temperature gradient in the basal ice q_c (Egholm et al., 2012a):

$$\begin{array}{}\text{(9)}& {\dot{m}}_{\mathrm{b}}={\displaystyle \frac{q_{\mathrm{b}}-q_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{ice}}{L}_{\mathrm{i}}}},\end{array}$$

where L_i is the latent heat for fusion of ice and ρ_ice is the density of ice (Table 1).

Figure 2 Paleo-topographic and bathymetric reconstructions. Panel (a) shows the PREQ experiment and (b) the MLQ experiment. Panels (c) and (d) show the differences between their respective panels above and the modern-day topography and bathymetry.

2.1.2 Topography and bathymetry

The focus of this study is to examine the influence of bed topography on ice-sheet behavior, exemplified by simulating the SIS on landscape configurations representing different periods in the Quaternary. For comparison, we simulate the SIS on modern-day topography and bathymetry over the last glacial cycle in a reference model. The reference experiment uses the global DEM GEBCO 2022 grid (GEBCO Bathymetric Compilation Group, 2022) sampled at 10 km × 10 km for the ice model. (The same grid resolution is used in all experiments.) Because of computational limitations, a model resolution higher than 10 km is not feasible. Having a higher resolution would allow us to resolve glacial morphology in higher detail and could lead to interesting findings regarding the influence of fjord systems in western Norway on ice-sheet dynamics. Here, we focus on larger features, such as the Norwegian Channel, where a 10 km resolution is sufficient. Throughout the model simulations, ice-driven isostasy is handled with a 2D uniform thin elastic plate model (e.g., Pedersen et al., 2014).

The pre-glacial landscape is adopted from Pedersen et al. (2021) and reconstructed using a source-to-sink approach that also considers (i) a component of glacial erosion that has taken place on the inner shelf, (ii) erosion-driven isostasy, and (iii) a component of dynamic topography (Pedersen et al., 2016). For further details on the approach see Pedersen et al. (2021). Here, we extend these previous reconstructions and remove the Quaternary sediment package from all sectors of the North Sea to reconstruct a realistic pre-glacial bathymetry for the entire region (Binzer et al., 1994; Rise et al., 2005; Nielsen et al., 2008; Knox et al., 2010; Gołędowski et al., 2012; Lamb et al., 2018). These additional sediment volumes, from outside of the Norwegian and Danish sectors, are not included in the landscape reconstruction at onshore Scandinavia. The result is a landscape representing a pre-glacial state before any major glaciations in Scandinavia, featuring a large submarine depression in the North Sea and a much narrower continental shelf along the Norwegian margin than at present (Fig. 2a, c). In addition to the PREQ experiment, two sub-experiments are presented: “PREQ-onshore” where only the onshore fjord erosion has been reconstructed (material added to present-day topography) and “PREQ-offshore” where only the offshore deposition has been reconstructed (material removed from present-day bathymetry). Neither of these additional sub-experiments considers the offshore sediment wedge on the shelf. With the sub-experiments we can assess which processes control the behaviors and ice-volume changes observed in the PREQ experiment.

For the middle–late Quaternary experiment, we reconstruct the bathymetry by estimating the volumes of erosion that have been carved into the modern-day seabed by ice streams on the Norwegian shelf and in the Norwegian Channel (Fig. 1). This bathymetric erosion is estimated using the geophysical relief method (e.g., Steer et al., 2012; Pedersen et al., 2021) on the present-day GEBCO 2022 global DEM (GEBCO Bathymetric Compilation Group, 2022), using a grid resolution of 1 km × 1 km and a sliding window radius of 35 km. The resulting filled bathymetry, which also fills fjords to sea level, is adjusted with the flexural isostatic response to loading using gFlex 1.1.1 (Wickert, 2016) with an effective elastic thickness of 15 km. This reconstruction of the Scandinavian morphology is meant to represent a state before the formation of the Norwegian Channel (Fig. 2b, d) and could represent an age of approximately ∼0.5 Ma. This approximate age is supported by the presence of buried mega-scale glacial lineations and drumlins in stratigraphic sequences of the North Sea suggesting that grounded ice has been present since ∼0.5 Ma, whereas the lack of these features in the older strata indicate that early Quaternary glaciations did not ground but only supplied icebergs to the North Sea (Dowdeswell and Ottesen, 2013; Rea et al., 2018).

3 Results

In this section we start by presenting the results from our reference model simulating the evolution of the SIS on the present-day topography and bathymetry over the last glacial period. Then we present the results of our two experiments with reconstructed topography and bathymetry, as well as how they differ from the reference model. Lastly, we present our findings regarding the formation of the Norwegian Channel.

Figure 3 Model output from three time slices of the reference experiment. Left column: early glaciation (72 ka); middle column: late–intermediate glaciation (22 ka); right column: glacial maximum (17 ka). Panels (a–c) show ice thickness, (d–f) mass balance, (g–i) depth-averaged deformation velocity, and (j–l) sliding velocity.

3.1 Reference model

To illustrate the spatial and temporal development of the SIS in our model simulations, we present model output from three snapshots in time (Fig. 3): minor ice build-up during early glaciation (72 ka), moderate glacial build-up during intermediate times of the glaciation (22 ka), and the glacial maximum that happens in these simulations at 17 ka. We note that the delayed timing of the glacial maximum in our models compared with the timing of the reconstructed maximum extent in Scandinavia (∼21–19 ka; Hughes et al., 2016) is a direct consequence of the chosen climate forcing, utilizing a glaciation index that peaks at 18 ka. We do not intend here to match the exact timing of the maximum extent (i.e., the LGM). During our simulated early glaciation, ice extent is limited to mountain regions with high topography and high latitude regions in Norway and Sweden (Fig. 3a). Mass balance is positive ∼1.5 m yr⁻¹ in high altitude regions at the Norwegian coast where precipitation is high and temperatures are low (Fig. 3d). Ice deformation and sliding are high (up to >50 and >200 m yr⁻¹, respectively) during early glaciation (Fig. 3g, j) in areas where the thin ice cover is controlled by the underlying topography which includes mountainous regions dissected by fjords and valleys.

During the intermediate glaciation, the ice sheet has advanced onto the shelf region, with grounded ice on the Norwegian margin, and the ice sheet has started to advance into the North Sea through the inner part of the Norwegian Channel (Fig. 3b). The mass balance reaches ∼1 m yr⁻¹ at the west coast of Norway (Fig. 3e), with values <0.5 m yr⁻¹ across most of the ice sheet and negative mass balance at the southwestern margin reaching $\sim -\mathrm{2}$ m yr⁻¹ where the ice is thin and velocities exceed ∼200 m yr⁻¹ (Fig. 3h, k). Along the coastal margin to the west, the mass balance is negative in a narrow zone where floating ice is melting fast. Sliding is notably high, reaching >200 m yr⁻¹ in the inner parts of the Norwegian Channel (Fig. 3k). The ice flow is still steered by topography in the high regions of southern Norway and in the Bothnic Bay, whereas the main divide in northern Scandinavia has shifted east, being largely independent of the underlying topography (Fig. 3e).

During the glacial maximum, the ice sheet reaches a thickness of >3000 m in the central parts (Fig. 3c) with a relatively low positive mass balance along the west coast of Norway (<1 m yr⁻¹; Fig. 3f) with the same general spatial pattern in accumulation and ablation as the intermediate glaciation (Fig. 3e) across the ice sheet. Sliding is high along the northwestern margin of the ice sheet (>200 m yr⁻¹) especially near the shelf break where ice is funneled towards the deeper ocean (Fig. 3l). For a while (∼5000 years) during the maximum expansion, the ice sheet merges with the BIIS in the western part of the North Sea, simulated as an ice wall (Fig. 3f, l). At that time, the ice flow rearranges into a divergent pattern from the ice saddle that emerges between the BIIS and the SIS. Consequently, the ice flows across the Norwegian Channel during the maximum extent instead of being focused in the channel itself, as the ice is diverged southward, driven by the surface slope of the ice sheet under this ice configuration (Fig. 3l). It is worth noting that the reference model captures a realistic placement of the LGM ice divide (Fig. 3f) in accordance with geological observations (Fig. 1; Olsen et al., 2013). Additionally, the ice divide of the saddle across the North Sea during the glacial maximum, when the SIS merges with the BIIS, closely resembles the ice divide suggested by Clark et al. (2022) using a combination of observations and modeling techniques. The glacial maximum ice extent in our reference experiment is within the maximum LGM ice extent (Fig. 1; Hughes et al., 2016), albeit with less ice towards the southern margin and more ice in the northeast.

Figure 4 Advance and retreat of the SIS in the reference experiment: (a) ice advance in 5 kyr intervals between model years 72 and 17 ka, and (b) ice retreat in 1 kyr intervals from 17 to 12 ka.

Build-up of the SIS from early mountain glaciation to the glacial maximum happens gradually with grounded ice on the Norwegian shelf forming 10 000 model years before the glacial maximum, and ice advance in the North Sea occurs over just 5000 model years approaching the glacial maximum extent (Fig. 4a). In contrast, the ice retreat is rapid with ice mass loss from the glacial maximum back to a state similar to that of early glaciation happening over just 5000 model years (Fig. 4b).

Figure 5 Differences in ice thickness for the (a–c) PREQ and (d–f) MLQ experiments compared with the reference experiment. Blue colors mean more ice in this experiment than the reference experiment and red colors mean less ice.

Figure 6 Panel (a) shows ice volume over time for the different experiments and (b) shows volume differences between the different experiments and the reference experiment. The dashed black line in (a) is the reference model, the red line is the middle–late Quaternary experiment, and the yellow line is the early Quaternary experiment. The green and blue lines represent the two sub-experiments of the early Quaternary experiment (offshore and onshore landscape changes compared with present day, respectively).

3.2 Results from PREQ and MLQ

In the model simulation representing ice-sheet behavior on an early Quaternary landscape morphology (PREQ; Figs. 2a, 5a, b, c), the ice sheet initially extends further than the reference model (Fig. 5a, purple color), particularly towards the Norwegian coast. At the intermediate stage (Fig. 5b), the ice sheet shows a smaller extent and thickness towards the Norwegian margin (Fig. 5b, orange color), whereas the ice extends further towards the south (Fig. 5b, purple color) with an ice thickness increase of >500 m in some regions. The location of the present-day Norwegian Channel shows a much thinner ice since this bathymetric depression is not present in the PREQ landscape reconstruction (Fig. 5b). At the maximum extent, the ice sheet is smaller both along the western and the southwestern margins (Fig. 5c, orange color), with a general decrease in ice-sheet thickness compared with the reference model (Fig. 5c, red colors). The reduced extent and ice thickness during the maximum extent result in ∼10 % lower maximum ice volume than the reference model (Fig. 6, orange curve). The large difference in ice volume between the PREQ experiment and the reference experiment is largely driven by differences in bathymetry (PREQ-offshore; Fig. 6a, green curve) as changes in topography do not lead to significant differences in ice volume compared with the reference model (PREQ-onshore; Fig. 6a, blue curve).

For the MLQ simulation that represents ice flow on a landscape morphology that existed prior to extensive erosion of the bathymetry by ice streaming (Figs. 2b, 5d, e, f), the ice sheet also starts slightly larger (Fig. 5d, purple color) compared with the reference model. At the intermediate stage, the ice sheet has already extended all the way across the North Sea (Fig. 5e, purple color), showing also a significantly thicker ice sheet in the adjoining regions at onshore Scandinavia. This trend is continued during the maximum extent, where the MLQ ice sheet extends even further, particularly towards the south (Fig. 5f, purple color). In general, the extent of the MLQ ice sheet is not changed along the Norwegian margin, where the width of the shelf has not changed for this simulation (Fig. 5e, f). The increased ice extent and ice thickness in the MLQ simulation result in a maximum ice volume that is ∼25 % more than that of the reference model during the intermediate stage and ∼5 % during the glacial maximum as a direct result of the changed bathymetry (Fig. 6, red curve).

3.3 Sliding in the Norwegian Channel

The erosive power of ice is a product of ice flux over a region with grounded ice (Patton et al., 2022) and is strongly correlated with ice sliding velocity (Cook et al., 2020), which means sliding velocity can be considered a proxy for erosive potential. Here, we explore whether our higher-order ice-sheet model can capture the erosive potential through sliding in the Norwegian Channel in the present-day bathymetry of the reference model and whether the model can predict erosion when the channel is not there in the PREQ and MLQ experiments. The ice dynamics in our reference simulation show significant sliding in the Norwegian Channel in four distinct phases (Fig. 7a, b, c, d). In the early glacial stage, the ice is sliding fast southeast of southern Norway along the deepest part of the channel (Fig. 7a). As the ice approaches maximum extent, the sliding pattern changes because of the different ice flow patterns that arise as an ice saddle emerges in the North Sea when the SIS merges with the BIIS (Fig. 7b). At this stage, ice flows south across the channel from the southern mountains of Norway, following the steepest surface gradient of the ice sheet. Instead, sliding is now mostly concentrated in the outer parts of the Norwegian Channel close to the North Sea Fan (Fig. 7b). During retreat, ice sliding continues in the outer parts of the channel but also becomes prominent along the southern tip of Norway with ice sliding towards the southeast, and in the inner parts of the channel near Oslo Fjord (Fig. 7c). Finally, as the ice sheet retreats further, continued sliding towards the North Sea Fan is complemented by a phase of southward sliding in the channel along the southwestern coast, a region that had not seen significant prior sliding (Fig. 7d). In Fig. 7e–h we show the same time slices for the MLQ experiment. Here, the ice extends further towards the west and has already formed a saddle between the SIS and the BIIS during the initial phase of the glacial cycle (Fig. 7e), and sliding is high towards the shelf break in the region that will later become the outermost part of the Norwegian Channel. Sliding velocities towards the shelf break are consistently high throughout the model simulation (Fig. 7f, h), whereas sliding accelerates in the inner parts of what will become the Norwegian Channel during ice retreat (Fig. 7g). In the last time slice, sliding velocity is lower than in the reference experiment but has the same general pattern (Fig. 7h), with sliding in some regions along the west coast of southern Norway. In Fig. 7i–l we present the time slices for the PREQ experiment. Across all four panels the patterns differ from the reference and MLQ experiments. Instead, we observe high sliding velocities towards the west across where the channel is today (Fig. 7i, j, k). In the last time slice we observe very little sliding as the ice has retreated mostly onshore at this time in the PREQ experiment (Fig. 7l).

Figure 7 Sliding velocity in southwestern Norway for reference model year (a) 23 ka, (b) 17 ka, (c) 14 ka, and (d) 13 ka. Same for the MLQ experiment (e–h) and the PREQ experiment (i–l).

4 Discussion

4.1 Ice extent and volume

The ice volume in our reference experiment reaches 6.5 M km³ at the glacial maximum, which is within estimates of SIS and Eurasian ice-sheet volume from previous studies (e.g., Hughes et al., 2016; Patton et al., 2016; Simms et al. 2019). The ice divide of the SIS in the reference experiment is in good agreement with observations (Figs. 1, 3f) which also affirms that our model captures an adequate representation of the ice sheet during the last glacial period. The differences in maximum ice extent between our reference experiment and observations (see Fig. 1; Hughes et al., 2016) can be attributed to the simple mass balance implemented in our model, using linear gradients, that does not capture the complex nature of the regional climate during the last glacial cycle but is an adequate approximation for our purposes. Geological observations suggest that the main ice advance in Denmark approaching the glacial maximum between 20 and 22 ka came from the northeast bringing till deposits of Middle Swedish provenance (Houmark-Nielsen, 2004), whereas the main ice advance into Denmark in our reference experiment comes from the north (Fig. 4a). A possible reason that our model does not capture this dynamic in the southerly ice advance could be the lack of subglacial hydrology in the model which can increase sliding rates (Egholm et al., 2012a). It could also be the lack of a more complex stress-dependent ice viscosity, where the Glen's flow law stress exponent can increase to n≈4 in some areas, which can increase the flow velocity by an order of magnitude (Millstein et al., 2022). These effects could be important especially in the southern parts of the ice sheet where the ice is thin and fast flowing during advance (Fig. 3c, i). Here, an even faster and thinner ice might be more sensitive to the low-relief topography of southern Scandinavia leading to a more westerly ice flow from Sweden into Denmark in agreement with the observations.

We cannot directly compare the ice extents in our experiments with reconstructions of past SIS extent as we use the same climate forcing between experiments, but we can assess whether differences in past ice-sheet extents follow the same trends as we see in this study that are based solely on differences in morphology. Batchelor et al. (2019) use empirical data to evaluate past Northern Hemisphere glacial extents and suggest best-estimate maximum southern extents of the MIS 12 (429–477 ka), MIS 16 (622–677 ka), and MIS 20–24 (790–928 ka) ice sheets to be somewhere between the best-estimate maximum MIS 6 extent and the LGM extent (132–190 ka; Fig. 1, dashed red line and black line), although the MIS 16 and MIS 20–24 maximum ice-sheet extents are highly uncertain. These reconstructions are based on very limited observations and, in some cases (e.g., MIS 12 and MIS 16), the estimates are mostly based on similarities in the δ¹⁸O curve (Batchelor et al., 2019). We show with this study that purely morphological differences in bathymetry between the last glacial period and ∼0.5 Ma (MLQ experiment, similar in time to MIS 12 and MIS 16) allow for larger ice-sheet extents simply owing to geomorphic changes during this time period. This suggests that both climatic and topographic forcings might have caused these (possibly) large ice extents of the middle–late Quaternary (MIS 12, MIS 16, and MIS 20–24). Indeed, our results showcase that a smooth bathymetry in the North Sea region (i.e., lacking glacial morphology), such as before the inception of the Norwegian Channel, could lead to earlier and more extensive southerly ice advance within a glacial period (Fig. 5e). On the other hand, our simulation of early Quaternary glaciations suggests that ice build-up across the North Sea was not plausible at this early stage of glacial landscape evolution. Indeed, in the PREQ experiment we found that the SIS could extend no further than the continental shelf during the early Quaternary (Fig. 5b, c). This is consistent with a study of buried glacial landforms in the central North Sea documenting iceberg plough marks in early Quaternary sediments (Dowdeswell et al., 2013; Rea et al., 2018). Our reconstructed early Quaternary ice sheet would have supplied icebergs that created these plough marks.

The differences we find in ice volume at the maximum glacial extent (∼5 % higher for MLQ, ∼10 % lower for PREQ) illustrate how differences in morphology affect ice volume independently of climate forcing. This has implications for the proxies we use for ice-volume history. Clearly the effect of glacial morphology explored here is local in nature, whereas the LR04 Benthic Stack we use as a glacial index and a proxy for ice volume is a global proxy. In addition, local ice volume also depends on global atmospheric circulation patterns which can lead to asynchronous development of the ice sheets during a glacial period (e.g., Liakka et al., 2016) that will also influence ice-sheet volume between glacial cycles. But landscape evolution has also played a significant role along other ice-sheet margins through the Quaternary, for example, leading to increased ice-sheet advance across marine sectors of the Antarctic ice sheet (Hochmuth and Gohl, 2019; Hochmuth et al., 2020). It should also be noted that the lack of ice shelves in our model could have a significant impact on grounded ice volume as buttressing effects of ice shelves can stabilize and advance grounding lines across the marine sectors of an ice sheet (e.g., Gasson et al., 2018). Nevertheless, according to this study, landscape morphology alone can account for up to ∼10 % difference in ice volume between glacial cycles for the Scandinavian region (∼25 % during ice build-up), implying that glacial landscape evolution could be an overlooked mechanism impacting local and global ice volume and thereby the interpretation of δ¹⁸O curves. This emphasizes the added uncertainty of landscape morphology on Quaternary ice-sheet reconstructions.

4.2 Formation of the Norwegian Channel

It is uncertain how and when the Norwegian Channel was formed, with studies estimating the time of formation to be between ∼0.35  and 1.1 Ma – with more recent studies suggesting younger ages (e.g., Sejrup et al., 2003; Hjelstuen et al., 2012; Løseth et al., 2022). In this study, we have assumed that the entirety of the Norwegian Channel formed after ∼0.5 Ma (MLQ). For the last glacial cycle, it has previously been proposed that the Norwegian Channel Ice Stream (NCIS) was active in stages but mainly during the LGM (e.g., Sejrup et al., 1998, 2003). According to an earlier study (Sejrup et al., 2016), ice streaming in the outer parts of the channel near the shelf break started close to the LGM with increased activity promoting ice retreat around 19 ka because of increased ice mass loss. The retreat translated southwards over time as the SIS unzipped from the adjacent BIIS after which ice streaming was mostly confined to the main trunk of the channel (Sejrup et al., 2016). A previous modeling study also suggests that the NCIS was active in stages with streaming in the inner parts of the channel leading up to, and became inactive during, the glacial maximum because of the saddle forming from the merging of the BIIS and the SIS (Boulton and Hagdorn, 2006). We find in our reference model with present-day bathymetry that ice streaming was active in the inner parts of the channel before the saddle formed between the BIIS and the SIS, after which ice streaming velocity increased dramatically in the outer parts of the Norwegian Channel near the shelf break and decreased in the inner parts of the Norwegian Channel as the saddle formed, consistent with other studies based on observations of, for example, subglacial landforms combined with dated sediment cores (Sejrup et al., 2016). On the other hand, our reference experiment does not mimic at any time an NCIS spanning the entire trunk of the Norwegian Channel, which would significantly contribute to ice mass loss from rapid grounding line retreat as is supported by observations (Sejrup et al., 2016). However, with this model setup we cannot rule out the occurrence of continuous ice streaming in the entire Norwegian Channel after the LGM. Indeed, some processes central to reproducing realistic ice stream behavior are not included in iSOSIA, such as enhanced basal melt owing to basal friction, leading to accelerated thinning in regions with rapid ice sliding as well as effects of internal friction and temperature advection on ice viscosity which can greatly amplify sliding velocities (Millstein et al., 2022; Bondzio et al., 2017). These mechanisms could contribute to highly elevated sliding velocities, especially in the NCIS, and could facilitate a propagation of the streaming activity we observe in the outer parts of the channel to the inner parts. In addition, the static ice wall we use to simulate the merging SIS and BIIS introduces a highly persistent ice saddle, which may introduce unrealistic streaming patterns and ice extent during NCIS retreat (Fig. 7c, d, g, h). Indeed, a previous study facilitates the retreat of the Norwegian Channel, with a negative surface mass balance anomaly in the southern sector of the North Sea, in order to match the ice margin to empirical reconstructions (Gandy et al., 2021).

Despite the Norwegian Channel being filled with sediment in the reconstructed bathymetry of our MLQ experiment, we find an ice streaming pattern that is comparable to that of the reference model for several parts of the model (Fig. 7a, b, c, d, e, f, g, h). Specifically, in the MLQ experiment, high sliding velocities are also present in what will become the inner part of the Norwegian Channel as the ice begins to advance offshore (Video 3 in the Video Supplement), although the ice velocity will be less focused compared with the reference model where the depression of the Norwegian Channel steers the ice even further (Fig. 7a). We stress, however, that because the ice advances faster offshore in the MLQ experiment, this sliding in the inner parts of what will become the Norwegian Channel happens prior to 23 ka (Fig. 7e; Video 3 in the Video Supplement). The MLQ experiment also shows high sliding rates where the outer part of the Norwegian Channel will form towards the shelf break (Fig. 7e, f, g, h), even extending further back in time than the reference experiment (Fig. 7a, e). This steering of ice towards the north–northwest in the MLQ experiment, steering which takes place before a bathymetric depression is formed, is mainly controlled by the steeper ice-surface gradient that occurs towards the shelf break in this simulation, when the ice advances into the offshore region and approaches the shelf break much earlier than in the reference experiment. This ice flow pattern begins before the saddle between the BIIS and the SIS formed but is amplified further by the ice saddle that forms in the North Sea as the ice cannot advance further towards the west (Video 3 in the Video Supplement). Our models can thus explain the initial formation of the Norwegian Channel in the innermost and outermost parts, starting from a bathymetry that had no prior imprint of the present-day channel. The MLQ experiment also shows sliding in other parts of what will become the Norwegian Channel later in the model simulation (e.g., Fig. 7g, h). However, we find these results less robust owing to the limitations of our model setup during the deglaciation.

On the other hand, the PREQ experiment showed no ice flow or sliding patterns similar to those of the reference model in the region that would later become the Norwegian Channel. Indeed, ice flow and sliding are at all times perpendicular to the future Norwegian Channel because of the sediment wedge that existed along the Norwegian coast and a steep ice-surface gradient towards the North Sea, sustained by the deep bathymetry of the North Sea that prevented grounded ice. Therefore, we find it likely that the carving of the Norwegian Channel could not have been initiated before the North Sea basin had been sufficiently filled with sediments. Instead, we find it plausible that the Norwegian Channel formed during multiple glacial periods since ∼0.5 Ma, consistent with a recent study indicating that the channel was formed prior to ∼0.35 Ma (Løseth et al., 2022). Our results are also in agreement with studies on the North Sea Fan (NCIS depocenter), suggesting that 90 % of the sediments in this fan are younger than ∼0.5 Ma (Hjelstuen et al., 2012).

5 Conclusion

We used a higher-order ice-sheet model to investigate the effects of landscape morphology on the SIS evolution and dynamics. Three different experiments where conducted: (i) a reference experiment resembling the last glacial cycle using modern-day topography and bathymetry, (ii) a middle–late Quaternary experiment with glacial morphological features in the present-day bathymetry filled with sediment, and (iii) a pre-Quaternary experiment simulating the SIS on a reconstructed pre-glacial topography and bathymetry. We found in the MLQ experiment that removing glacial morphological features in the bathymetry allowed for faster and further southward expansion at similar climatic conditions allowing for a larger ice sheet. Contrary to this, we found in the PREQ experiment that the early Quaternary bathymetry did not allow for the SIS to advance as far westward and southward, thereby limiting the size of early glaciations and preventing a merging of the BIIS and the SIS. Looking at the prominent glacio-morphological feature, the Norwegian Channel, we found that the PREQ experiment did not allow for significant ice streaming in this area and that the channel was more likely formed after the North Sea was filled in with glacial sediments. Furthermore, our results suggest that ice streaming occurred in distinct stages along the trunk of the channel with high ice sliding in the inner parts before the LGM and sliding in the outer parts of the channel close to the shelf break during the LGM. Our results also show that sliding in the inner parts of the channel ceased because of divergent ice flow when the BIIS and the SIS merged and formed a saddle across the North Sea.