Younger Dryas transgression in western Norway: a modelling approach

Fennoscandia has experienced major uplift since late-glacial time. However, in the Younger Dryas (~12,700–11,500 calibrated years before present) 10 metres of transgression occurred at the Norwegian west coast. This transgression was most significant in the area around Bergen; both to the north and to the south of this area the relative sea level fell during the Younger Dryas. Thus, in the area of transgression the shorelines are twisted; the 60 m Younger Dryas isobase crosses the 60 m Allerød isobase. The same section of the coast that was transgressed during the Younger Dryas also experienced a major ice-sheet readvance. In part of this area the ice sheet readvanced by at least 40 km. It seems likely that there is a causal connection between the ice readvance and the shoreline transgression. In this paper we report modelling of the isostatic and eustatic response of a readvancing Younger Dryas ice sheet. It is found that a halt in the glacial isostatic rebound is the major cause for the observed transgression. This, however, presupposes an Earth rheology with a low-viscosity asthenosphere and a weak effective elastic lithosphere.


Introduction
The Earth has been and is still readjusting after being subjected to huge ice loads during the last glaciation. The general picture in Norway is that the sea level has been steadily falling through late-glacial time. However, a significant relative sea-level rise during the Younger Dryas (YD) constitutes a conspicuous feature of the Late Weichselian sea-level history of western Norway. One of the first documentations of the late-glacial transgression in western Norway was from the island of Sotra, west of Bergen (location; cf., Fig. 1B). The early reconstructed sea-level curve of Sotra was based on data from a large number of isolation basins Krzywinski & Stabell, 1984;Anundsen, 1985).
The so-called isolation basin method (Hafsten, 1960) is based on age dating of stratigraphic boundaries between marine and lacustrine sediments in lakes corresponding with the isolation of the lake from the sea, which is the time when the sea level fell below the level of the outlet threshold of the lake basin. The position of the sea level relative to the outlet of the basins determines whether the sediments are marine, brackish or lacustrine. The age of the sea level corresponding to the outlet of each lake is determined by identifying and dating the boundary between the brackish and the lacustrine sediments in the cores. A relative sea-level curve with possible transgressions and regressions can be constructed with data from a number of isolation lakes of various altitudes within a limited area (see more details in Kjemperud, 1981;Svendsen & Mangerud, 1987;Fjeldskaar & Amantov, 2017). Since the first documentation of the Sotra curve, the curve has been significantly improved by collecting new cores and better analysis (Lohne et al., 2007(Lohne et al., ). et al., 2007Lohne et al., 2007). Numerical modelling carried out in the 1980s suggested that the transgression was caused by the YD glacial readvance in western Norway (Mangerud, 1977(Mangerud, , 2004Fjeldskaar & Kanestrøm, 1981;Anundsen & Fjeldskaar, 1983). The numerical modelling indicated that the glacial regrowth slowed down the glacial isostatic rebound and that the growing ice-mass caused a geoid effect (Fjeldskaar & Kanestrøm, 1980;Anundsen & Fjeldskaar, 1983). Lohne et al. (2007) found that the YD shoreline had a similar gradient or was slightly steeper than the Allerød shoreline, which indicates that the isostatic rebound ceased (or possibly was slightly reversed) during this time span.
We have here revisited the YD transgression of Sotra as well as the regional YD high-stand shorelines by new numerical modelling. The modelling is based on a high-resolution deglaciation model combined with a high-resolution glacial isostatic model characterised by a thin, effective elastic lithosphere and a low-viscosity asthenosphere.

Observations
The Norwegian west coast has experienced a significant sea-level transgression in the Younger Dryas (YD; ~12,700-11,500 BP) from the Allerød (AL; ~13,900-12,700 BP) low-stand (Fig. 1A). The transgression was most significant in the outermost parts of western Norway in the areas between Sognefjord and Boknafjord, where the transgression from AL to YD was close to 10 m. Both north of the Sognefjord area and south of the Boknafjord area (for locations, see Fig. 1B) the relative sea level fell during YD. In the area of transgression of the YD and AL shorelines the 60 m YD isobase crosses the 60 m AL isobase (Fig. 1B). The AL shoreline shows increasing height towards the inland, with a gradient of 1.3 m/km; and the gradient seems to have been stable from AL to YD (Lohne et al., 2007). The reconstructed relative sea-level curve for Sotra ( Fig. 1A) is based on emerged lake basins covering the last 14,500 years. The stratigraphy and the age datings from three basins (Gardatjønn, Sekkingstadtjønn and Hamravatn) define the YD transgression. The AL lowstand is determined at 29.6 m a.s.l., which seems to have been more or less stable at this site for more than 500 years. After this still-stand, the relative sea level started to rise and reached the peak of the YD transgression at 39.3 m a.s.l. at some time prior to 11,500 BP. The magnitude of the YD transgression was 9.7 m at Sotra and the average sea-level rise during the transgression was 7 mm/ yr (Lohne et al., 2007).
A YD transgression has been determined at several locations on the southwestern coast of Norway (Anundsen, 1985;Helle, 2004;Lohne et al., 2004;Helle et al., 2007). No positive evidence for significant lateand post-glacial faults has been found in the Sotra region (Lohne et al., 2007). The regional extent of the transgression points to a mechanism of a different nature than fault movements.
The same section of the coast that was transgressed during the YD also experienced a major YD ice-sheet readvance (Fig. 1B). In part of this area the ice sheet readvanced by at least 40 km . Mangerud et al. (2016) have presented 90 radiocarbon dates from 36 sites showing that a glacial readvance took place in western Norway during the YD. The late culmination of the advance is accurately dated to the end of YD (~11,500 BP). A time-distance diagram for the Hardangerfjord area (for location, see   Lohne et al., 2007). Larsen et al. (2010). It is worth mentioning that this evolution is compatible with the reconstructions of the recently published time-slice reconstruction by . The AA1 ice model readvance from AL to YD for western Norway is based mainly on Lohne et al. (2007).
To our knowledge, AA1 is the only deglaciation model with high enough resolution to properly model the effect of a YD glacier regrowth. Neither of the commonly publicly available ice models ICE-5G (Peltier, 2004) or ICE-6G (Argus et al., 2014;Peltier et al., 2015) is of high enough resolution, and both show a continued decay of ice volumes in western Norway from the AL to the YD. These deglaciation models are therefore not relevant as input for modelling the effect of the YD glacier regrowth.
The method used to compute ice-sheet thickness development of the AA1 consists of the following steps (more details are given in the Discussion section below and in Amantov & Fjeldskaar, 2013):  Estimation of a general ice sheet with averaged values and shapes associated with a viscoplastic flow of ice known from present-day ice sheets,  Modification of ice-thickness distribution and different ice-flow velocities due to subglacial topography,  Corrections of ice thickness by ice streams, topographic roughness at the ice base, areas of discharge, etc.
The full reconstructed deglaciation history of the AA1 used in the calculations is shown in Fig. 3. The decay (or

Modelling
The Earth's response to glaciers is here modelled by using a flat Earth layered viscous model overlain by an elastic lithosphere. For more details on the modelling technique, see Fjeldskaar (1994Fjeldskaar ( , 1997. Fjeldskaar (1997) showed that the calculated isostatic equilibrium deflections by our flat Earth filtering technique show insignificant deviations from the deflections calculated for an isotropically elastic, uniformly thin, spherical shell for a load up to (at least) 1000 km in radius.
There are advantages and disadvantages with a regional flat Earth model. The main challenge for such a model is that it cannot provide detailed sea-level changes as a global model can. The advantage, however, is that the modelling can be done with high resolution and short CPU time. The spatial resolution in the modelling reported here is 10 km. This high resolution is a requirement for modelling of the regrowth of the YD glacier.

Deglaciation history
A realistic model of the deglaciation history is essential for glacial isostatic modelling. A model of the deglaciation history consists of two main elements: 1) a model of the glacial extent over time, and 2) a model of the ice thickness over time. We have used the deglaciation model AA1 from Amantov & Fjeldskaar (2013), which is developed independently from calibration of the glacial isostatic modelling. The spatial evolution of the ice sheet of AA1 from 20,000 BP (LGM) to 14,000 is based mainly Fjeldskaar & Amantov (2017) studied the YD (12,000 BP) shoreline diagrams of coastal Norway and concluded that an Earth model with a thin elastic lithosphere and a lowviscosity asthenosphere matches the observed shoreline gradient. The glacial isostasy in the present study is calculated with an Earth rheology model consistent with the above study; an effective elastic thickness of ~30 km (flexural rigidity 5 x 10 23 Nm), and an asthenosphere with a thickness of 75 km and a viscosity of 1.3 x 10 19 Pa s above a mantle of viscosity 10 21 Pa s. This Earth rheology is also strongly suggested based on previous studies of the Fennoscandian present rate of uplift (Fjeldskaar & Cathles, 1991). The method for calculating the isostatic response is described in Fjeldskaar (1994Fjeldskaar ( , 1997. The full deglaciation history (shown in Fig. 3) from the growth) from one ice-sheet configuration to the next is assumed to have uniform speed. An enlargement of the reconstructed ice thickness in western Norway for the AL and YD glaciers is shown in Fig. 4.

Glacial Isostasy
Observations of relative sea level (RSL) can be visualised in two ways: • Shoreline displacement curves, showing the vertical displacement at a certain location, • Shoreline diagrams, showing the displacement and tilting of paleo shorelines. LGM (Last Glacial Maximum) was used to calculate the glacial isostasy, but we here only present the resulting glacial isostasy from 13,000 BP (AL) to 11,600 BP (YD). The AL ice sheet was assumed stable over 600 years (from 13,500 BP to 12,900 BP) before the ice sheet readvanced over 700 years to the YD position at 12,200 BP. The YD ice sheet was assumed stable until it underwent fast melting which started at 11,600 BP. See Fig. 2 for a timedistance diagram that was used in the calculations.
Hydro-isostasy depends on eustatic (meltwater additions to the ocean) changes, but also on the paleo relief, because the extent of the area covered by water has not been constant through time. The paleo relief is again a function of the glacial isostatic response. The area covered with water (ocean or lake) depends on the extension of the glaciers. The spatial distribution of ocean water is combined with eustatic sea level in the calculations of hydro-isostasy. The global eustasy, as measured on the Barbados curve (Fairbanks, 1989), was 20-25 m in the period 13,000-11,600 BP (see discussion about eustasy below). These eustatic sea-level changes were applied to the spatial distribution of oceans (as shown in Fjeldskaar & Amantov, 2017, Fig. 5).
The resulting glacial and hydro-isostasy from 13,000 BP to 11,600 BP is shown in Fig. 5. The isostatic rebound significantly halts in parts of western Norway, especially between the Hardangerfjord and Sognefjord areas due to the glacial regrowth.

Elastic effect
In addition to the dynamic isostatic process, there might be an elastic effect of the ice growth. Due to the ice growth there will be an immediate elastic deformation, proportional to the stress when a force is applied to the Earth's surface. Almost all solid rocks behave elastically when the applied forces are not too large, and return to their original shape when the force is removed. A possible elastic effect of the additional load due to ice growth could add to the pattern of Fig. 5, because elastic deformation happens instantaneously by changes in loads. Fig. 6 shows the calculated elastic effect from 13,000 BP to 11,600 BP, using the method in Fjeldskaar (2000). Fig. 6 shows that the elastic effect of readvance in the area is less than one metre; elastic deformation will thus not contribute significantly to the sea-level transgression of this area.

Geoidal eustasy
Another possible effect of the regrowth is a potential geoidal rise due to an increased gravitational effect of the ice. Geoidal eustasy is defined as changes in the ocean water distribution caused by variation in the Earth's gravity field. Geoidal eustasy is an important eustatic factor, and will result in non-uniform sea-level changes over the globe. The three types of eustasy (Mörner, 1976) are: 1) glacial eustasy (caused by variation in the ocean water volume); 2) tectono-eustasy (controlled by variation in the ocean basin volume) and 3) geoidal eustasy. Factors 1 and 2 are giving global uniform eustatic changes, but due to hydro-isostasy the sealevel history will differ between oceanic islands and continental margins. An island moving with the sea floor will record the full sea-level change and is commonly approximated by the curve from Barbados. The Barbados curve (Fairbanks, 1989) has been used in this study as a measure of eustatic factors 1) and 2).
It has previously been demonstrated that the geoidal eustatic effect can be significant due to ice growth from AL to YD (Fjeldskaar & Kanestrøm, 1980). Our new modelling gives a geoidal effect of the glacial regrowth from 13,000 BP to 11,600 BP as shown in Fig. 7. This geoidal change in western Norway is due to two effects -glacial regrowth from AL to YD and movement of the solid Earth over the same period. When there is a slow   deglaciation (or regrowth) the related geoidal eustatic deflections tend to be smeared out, because lost (or increased) ice mass will be compensated by changes in crustal mass. In our area the geoidal effect could be ~3 m. This effect will be added to the hydro-and glacial isostasy given in Fig. 5 and the global eustatic effect (based on the Barbados curve). The method for the geoidal calculations is given in Fjeldskaar (1991).

Modelled total effect of glacial readvance
The calculated added effect of the glacial readvance, glacial and hydro isostasy combined with eustatic change, is shown in Fig. 8. The area of predicted YD transgression in western Norway is clearly noticeable and covers large parts of the region. The area of predicted transgression exceeding 10 m spans the area from the Hardangerfjord to slightly north of the Sognefjord area. The observed YD shorelines in the Sotra area show a tilt of 1.3 m/km, and based on Fig. 5 it is evident that there will be very small changes in the gradient from AL to YD; as observations also show (Lohne et al., 2007).
The predicted area of transgression is the same area that will experience a regional YD high-stand. As shown in Fig. 9 the modelling results are in good agreement with the observations.

Discussion
The uncertainties in our calculations are related to five factors: 1) GIA model; 2) neotectonics; 3) data on sea level and glacial history; 4) ice model; 5) eustasy.

The GIA model
In the present modelling we have used an Earth model with a uniform mantle of viscosity 10 21 Pa s overlain by a 75 km-thick low-viscosity asthenosphere of viscosity 1.3 x 10 19 Pa s and capped by a thin lithosphere with an effective elastic thickness of ~30 km (flexural rigidity 5 x 10 23 Nm).
There is, however, no agreement in the literature as to which rheology model best fits the data of postglacial rebound in Fennoscandia. Steffen & Wu (2010) give a review of modelling results on Fennoscandian rheology based on observations of postglacial and ongoing uplift in the area. The suggested Earth rheology for Fennoscandia consists of an effective elastic lithosphere thickness (t e ) between 75 and 160 km, and a viscosity of the lower mantle up to 100 times higher than the upper mantle. This is a reflection of the fact that the majority of present-day modellers relies on the so-called channel  flow model for glacial isostatic calculations (see Cathles, 1975 andAmantov, 2017 for a thorough discussion).
However, analysis of tilted paleo shorelines in Norway shows that models with a thick elastic lithosphere (80-160 km) and a considerable viscosity difference between upper and lower mantle give a significant mismatch to the observations (Fjeldskaar, 1994(Fjeldskaar, , 1997Fjeldskaar & Amantov, 2017). A best fit with the observed tilted paleo shorelines is achieved with an Earth model with a low-viscosity asthenosphere and a weak lithosphere. However, we have in the following tested glacial isostasy from the AL to the YD for three different Earth rheology models (Fig. 10); all with a thick elastic lithosphere.
(i) Fig. 10A; a low-viscosity asthenosphere and 120 km-thick elastic lithosphere. With a total eustatic global change of ~22 m, there could be a minor sealevel transgression from the AL to the YD over the entire western Norway. A 10 m transgression on the Sotra island with this Earth rheology would require a 5 m higher global eustatic change. However, a higher eustatic change will increase the mismatch with the observations as it will give significant transgression in the southernmost and northernmost areas where a falling sea level is observed from AL to YD (Fig.  1B).
(ii) Fig. 10B; an upper mantle viscosity of 0.7 x 10 21 Pa s, lower mantle viscosity of 2 x 10 21 Pa s, and a lithosphere thickness t e ~140 km (as favoured by Kierulf et al., 2014). A global eustatic change of ~22 m would not give a transgression in western Norway. A transgression of ~10 m at the Sotra island would require an at least 15 m higher eustatic change during this period, apparently excluding this Earth rheology as a viable option.
(iii) Fig. 10C; upper mantle viscosity 0.3 x 10 21 Pa s, lower mantle viscosity 5 x 10 21 Pa s, and a lithosphere thickness t e ~75 km (as favoured by Lambeck et al., 1998). With a eustatic change of ~22 m, this Earth rheology predicts a small transgression (~3 m) in western Norway. However, a 10 m transgression from AL to YD would require at least a 7 m higher eustatic change during this period. This Earth rheology model also predicts a significant transgression far north of the Sognefjord area, which does not match with the observations.
Two of the Earth rheology models could give a small YD transgression in western Norway. However, in addition to the mismatch in the southernmost and northernmost parts of western Norway, there are two additional important observations that seem to exclude these Earth rheologies; 1) the observed gradients of the YD shorelines of 1.3 m/km; 2) observations show that no tilting took place between the mid AL and the late YD (Lohne et al., 2007). Fjeldskaar & Amantov (2017) showed that models with a thick elastic lithosphere cannot match the observed tilting of the YD shorelines. Even if the models of Fig. 10 do not give significantly different shoreline tilt between YD and AL, these models are not viable options due to the calculated low gradients of YD shorelines in the area. Thus, there seems to be an upper bound on the elastic lithosphere based on the effect of the YD glacial readvance. The upper bound of the elastic lithosphere thickness is significantly less than 75 km. More discussions on the modelling uncertainties are given in Fjeldskaar & Amantov (2017).

Neotectonics
The historical seismicity in Fennoscandia (including western Norway) is remarkably high for an interior plate region, with high seismicity concentrated in specific regional areas. Fjeldskaar et al. (2000) found indications of neotectonic movements during the postglacial uplift of Fennoscandia. It is therefore not unreasonable to assume that parts of the observed western Norway YD high sea-level stand could be due to faulting (Helle et al., 2007). However, no positive evidence for late-or post-glacial faults has been found in the Sotra region, and the linear shape of the YD shoreline indicates that potential vertical fault movements are of minor importance (Lohne et al., 2007). There might, however, be local sea-level reconstructions in western Norway which carry components of neotectonic faulting (Helle et al., 2007), but in the present study we have assumed that the regional pattern of YD sea-level high-stand is not significantly influenced by fault movements.

Data on palaeo sea level and glacial history
The above-mentioned analysis of the YD shoreline diagrams (Fjeldskaar & Amantov, 2017) included shoreline data from the same location as our current area of interest, the Sotra/Bergen area. The gradient of the YD (12,000 BP) shoreline of Sotra is 1.3 (± 0.1) m/ km (Lohne et al., 2007), and the same (or slightly higher) for the AL. The gradients of the palaeo shorelines are constructed by dating deposits in isolation basins. This is the most precise method available to determine both the elevation and the age of former sea levels (see more details in Svendsen & Mangerud, 1987;Fjeldskaar & Amantov, 2017).
The sea-level curve of Sotra (including the YD transgression) is constructed with the same precise method. Both the reconstructed tilting and the determined transgression are assumed to be very precise, with limited uncertainty. A similar high precision is assumed to be the case for the reconstructed glacial history of western Norway. The regrowth of the ice in the YD is based on 90 radiocarbon dates from 36 sites. The uncertainty of the age-dating is less than 500 years (see more on the dating uncertainties in Mangerud et al., 2016 andLohne et al., 2007).  (Kierulf et al., 2014); (C) upper mantle viscosity 0.3 x 10 21 Pa s, lower mantle viscosity 5 x 10 21 Pa s, and a lithosphere thickness t e ~75 km (Lambeck et al., 1998). The predicted area of transgression is shown by the grey colours. The contour interval is 5 m. The red line shows the 10 m isoline.

The ice model
The outline of the glacier during the deglaciation is based on observations, but unfortunately the palaeo ice thicknesses are not possible to observe, with some limited exceptions (Mangerud et al., 2013). For the ice thickness we generally have to rely on modelling combined with calibration to current glaciers. Our deglaciation model is of high resolution, which is a requirement for calculating the effects of the YD glacial regrowth (cf., Figs. 3 & 4). Our estimations of the ice thickness consist of: a) preliminary initial assessment of an oversimplified general ice-sheet sketch with averaged typical values and forms known from current ice sheets, and in agreement with Glen-Nye's flow law. We are using approximate glacier mass-balance and separate volume control for the growing and decaying stages, and a variable long-term balance ratio between ablation and accumulation gradients. We have used outlines of the ice sheets with ~1000-year intervals for prediction of the ice thickness.
b) estimation of ice thickness distribution from a given subglacial topography. Further zonal corrections of the ice thicknesses are used with approximations of possible ice streams with variable stress at the base and small basal drag, and a variable substratum of the ice sheets over time and space due to sedimentation and erosion, areas of different termination, heat flow, etc.
The starting point of the modelling is a simplified type of general shallow-ice approximation model combined with a shallow-shelf approximation (Pollard & DeConto, 2009). This model includes parameters for basal resistance. Analysis of (and calibration of) the present ice thicknesses in Greenland and Antarctica was performed using the ETOPO1 global relief model of the Earth's surface and ice base compilation.
A mathematical fit of ice thickness variations was based on numerous 2D slices that were used for model calibration. The resulting simplified preliminary icesheet model undergoes further improvement and corrections as described in Amantov & Fjeldskaar (2013). Topographic ice streams are taken into account as regions with variable properties (Pattyn, 2003). We adjust zonal corrections according to variations in slipperiness at the base of the ice sheet, and also rate of accumulation -wastage balance between continental and oceanic segments -and in some cases slope gradient. Our procedure gave the AA1 deglaciation history (Fig.  3), including the AL and YD ice thicknesses (Fig. 4) and the reconstructed ice surface of the AL and YD as shown in Fig. 11. We have to bear in mind that there could be significant uncertainty in the modelled ice thicknesses, even with calibration to the Antarctica and Greenland ice sheets. It is, however, worth mentioning that the ICE-5G deglaciation model gave similar YD shoreline gradients as the ice model AA1 for an Earth model with thin elastic lithosphere and a low-viscosity asthenosphere (Fjeldskaar & Amantov, 2017). The ICE-5G (and ICE-6G) model, however, is not relevant for detailed modelling of the Younger Dryas regrowth since, as mentioned above, it does not take into account the more detailed observed glacial history of western Norway.
In opposition to the traditional view of a Younger Dryas readvance, several geoscientists have argued that the Hardangerfjord area was ice-free during the YD and that the deposition of the outermost moraine (Herdla-Halsnøy) must be older than YD (Helle et al., 1997(Helle et al., , 2007Helle, 2004). More recent results, however, have confirmed that the Herdla-Halsnøy Moraine is of YD age (see Lohne et al., 2012and Mangerud et al., 2013 for an in-depth discussion on this topic).

Eustasy
As mentioned above, we know that eustatic changes are not globally uniform (Fjeldskaar, 1989). It is generally assumed that the sea level at islands far from the coast will register the full eustatic change because the islands are moving isostatically with the ocean floor. At the coastal areas the observed sea-level change will be different, because the continents are uplifting in response to the subsidence of the ocean floor.
Here, we assume that the sea-level curve from Barbados registers the global post-glacial eustatic change. From 13,000 BP to 11,600 BP, a 20-25 m sea-level rise has been registered in Barbados. Geoidal eustasy, caused by gravity change, will redistribute the ocean water. We have above shown our calculation of the geoidal change in western Norway at the relevant period -due to movement of the solid Earth and regrowth of the ice (Fig. 7). Whether the calculated geoidal change due to glacial regrowth is the only geoidal eustatic effect that affected the sea level in western Norway from the AL to the YD is, however, uncertain.

Conclusions
Modelling shows that glacial isostatic response of a readvancing Younger Dryas ice sheet is the major mechanism for the observed transgression in western Norway. The glacial regrowth from the Allerød (AL) to the Younger Dryas (YD) led to a halt in the postglacial rebound, but also to a rise in sea level due to geoidal eustasy. Together with the general global eustasy, this could have caused a transgression at the island of Sotra of ~10 m. The transgression is of regional extent and has caused the 60 m YD isobase to cross the 60 m AL isobase in parts of western Norway.
This explanation for the observed sea-level transgression and related shoreline tilting in western Norway presupposes a regional Earth rheology with a low-viscosity asthenosphere, and an effective elastic lithosphere thickness of only ~30 km.