Introduction

The Laptev Sea and the East Siberian Sea represent large parts of the Russian Arctic coast, and cover the East Siberian Shelf. This continental shelf extends over the East Siberian Platform and the North American Plate and is characterised by a broad and shallow shelf sea, which extends hundreds km offshore with water depths of less than 100 m. The respective plate boundary cuts through the Laptev Sea in North–South direction (Fig. 1).

Fig. 1
figure 1

Setting of studied area (orange box) at Northern Hemisphere glaciation. Land surface at LGM (20 ka BP) is shown in green. Ice sheet thicknesses according to ICE-5G [27] are superimposed. The respective ice sheets are labelled in bold lettering. Present-day shoreline is shown in dark grey; plate boundaries according to Bird [6] are shown in black, separating here the Eurasian (left) and North American Plate (right); the rotation pole describing the relative motion between the East European platform and the North American plate according to Franke et al. [11] is shown by the black star south of the Laptev Sea

Full size image

The East Siberian Shelf is covered by vast areas of submarine permafrost reaching a thickness of up to several hundred metres (see [28]). Submarine permafrost as terrestrial relict mainly formed during glacial periods on arctic lowlands which were flooded by rising sea or by coastal retreat since then [26]. After inundation, subsea permafrost degrades from above due to thermal (warmer) and chemical (salt water infiltration) influence of sea water [24]. From below acts the geothermal heat flux as permafrost-degrading process [29]. Overduin et al. [25] propose that new permafrost can also form on very shallow shore faces in the presence of relict permafrost. The close relation between permafrost formation/degradation and local sea level demands a realistic reconstruction of the sea-level history when interpreting the detected submarine permafrost on the continental shelf.

Retracing the relative sea-level history of North East Siberia is hindered by the lack of precise sea-level indicators. The relatively large distance of the shelf to the Northern Hemisphere ice sheets (see Fig. 1) suggests to consider the eustatic sea-level change to reconstruct the late-glacial to Holocene transgression history. For instance, Bauch et al. [5] and Nicolsky et al. [22] applied the eustatic sea-level curves of Fairbanks [8] and Fleming et al. [10], respectively.

But, applying directly the eustatic sea-level change, the authors neglect regionally varying contributions due to glacial isostatic adjustment (GIA). The surrounding Northern Hemisphere ice sheets result in uplift during and after the termination in the formerly glaciated regions and subsidence in the surrounding forebulge regions. This process of glacio-isostasy is governed by the flexural behaviour of the lithosphere and the viscous flow of mantle material in response to ice loading. Furthermore, the process diminishes from the formerly glaciated regions and is often neglected in the far field.

Deformations due to the eustatic sea-level rise appear as hydro-isostasy around the continental margins [10, 20, 31]. The mechanism of hydro-isostasy can be described as the flexure of the lithosphere in response to the changing water load: The global sea-level drop of 120 m during the last glacial maximum (LGM) [10] had unloaded the ocean bottom which resulted in an isostatic uplift of the ocean basins by 40 m when measured against the continents which do not suffer from such a load change (we assumed here the density ratio between water and mantle material to be one-third). During the termination of the glacial ice sheets, the sea level rose again, resulting in a loading of the ocean basins. At the continental margin, the elastic flexure of the lithosphere controls the bending of the lithosphere due to this loading contrast. This flexure is, similar to post-glacial rebound, an ongoing process due to the retarded response of the viscoelastic mantle material with subsidence offshore and uplift onshore.

The hydro-isostasy at the continental shelf of East Siberia is further complicated as it was almost completely land fallen during the last glacial period [5] which shifts the flexure of the lithosphere towards the continental shelf margin.

Both mechanisms, hydro- and glacio-isostasy, generate local sea-level variations that systematically deviate from the eustatic signal and depend on the lithosphere and mantle structure.

Whereas the East Siberian Sea is located on the stable North American Plate, the tectonic setting of the Laptev Sea is more complicated. In the Eurasian Basin, the Gakkel Ridge represents the boundary between the North American Plate and the East Siberian Platform and propagates under the Laptev Sea, an assumption which is based on the interpretation of seismic profiles [7, 14]. Due to the location of the rotation pole of this rift zone some hundred km south of the Laptev Sea [11] (see star in Fig. 1), it is expected that the ridge is passive and the observed seismicity is driven by the extensional stresses of the surrounding tectonic plates and not by an active rifting [12]. This is confirmed by tectonic modelling (e.g. [1]) and the extremely low spreading rate in the eastern part of the ridge of less than 0.7 cm/year [19]. The markedly lower seismic attenuation along the rift axis on the continental shelf than in the Eurasian basin confirms the absence of hot material below the shelf [11]. Therefore, a weak lithosphere and low viscous asthenosphere below the Laptev Sea will be not considered in the present model set-up and a pronounced structural change from the East European Platform to the North American Plate is unlikely to exist.

In a conceptual way, we will discuss the impact of glacial isostatic adjustment in terms of the glacio- and hydro-isostasy on likely reconstructions of the transgression history in this region since LGM, which is based on an ensemble of possible viscoelastic earth structures. Partly, this forward-modelling study will help to identify possible locations where the reconstruction of the transgressive history by geological sampling is demanded.

Method

The deformational response of the solid earth to the changing ice and water load distribution is determined by the spherical symmetric viscoelastic earth model VILMA [18], where we assume an incompressible viscoelastic material law. The viscosity structure of the mantle is separated at the 670-km discontinuity into an upper and a lower mantle shell, and an elastic lithosphere of constant thickness is considered at the top. Viscosities of the lower and upper mantle are varied in a considerable range, whereas for the elastic lithosphere its thickness is varied. All models are considered as spherically symmetric. In total, we consider an ensemble of 70 different earth model parametrisations (see Table 1). The lateral resolution of the earth model is 120 km. This is sufficient for the considered loading scenario and for the considered lithosphere thicknesses. The lithosphere acts as a strong low-pass filter which reduces the viscoelastic displacements below 200 km wavelength markedly.

Table 1 Parameterisations for the considered earth model ensemble
Full size table

To keep the parameter space simple, we do not vary the glaciation history, but consider the last glacial cycle according to ICE-5G [27]. The load distribution for different time steps is provided on a global Gauss–Legendre grid with 256 × 512 grid points.

To consider the water and ice redistribution consistently, we apply the method of Hagedoorn et al. [13] to solve the sea-level equation introduced by Farrell and Clark [9]. In this solution method, sea-level variations are determined on the same grid as the glacial load. Therein, the change of the coastline due to the changing water level is considered as well as its loading effect. The formulation is comparable to that of Kendall et al. [16], where we iterate the present-day topography. In solving the sea-level equation, water mass conservation is considered as the gravitationally consistent redistribution of water in the ocean basin, which is deformed by the additional loading. This means, that hydro-isostasy is implicitly solved for when considering the sea-level equation. This fact might have hidden the importance of hydro-isostasy resulting in a deviation from eustatic sea level when considering relative sea-level change along continental coasts [17].

The solution of this calculus at each epoch, t, results in spectral representations of the vertical displacement field of the surface, \(U_n^m(t)\), and of the gravitational potential, \(E_n^m(t)\), as well as the constant shift of the potential surface to which the sea level will adjust, \(\Delta s (t).\)

From these fields, we determine the relative sea level at each location, \({ \Omega }= (\lambda , \phi )\), and epoch, t, according to

$$\begin{aligned} h_\mathrm{RSL}({ \Omega }, t)&= \sum _{n=1}^{n_{\max }} \sum _{m=-n}^{n} \left[ \, E_n^m (t) - U_n^m (t) \, \right] {Y_n^m}^{*} ({ \Omega })\nonumber \\&\quad + \Delta s (t) \, - \, h_\mathrm{RSL}({ \Omega }, t_\mathrm{pt}). \end{aligned}$$
(1)

In this study, the summation is performed over the spectral coefficients up to maximum degree, \(n_{\max } = 170.\) \({Y_n^m}^{*}\) represents the fully normalised spherical harmonic function of degree, n, and order, m (e.g. [30]). The palaeo topography follows directly:

$$ h_\mathrm{topo}(\Omega, t) = h_\mathrm{topo}(\Omega, t) \, - \, h_\mathrm{RSL}(\Omega, t) $$
(2)

The match of the left- and right-hand side at \(t=t_\mathrm{pt}\) is guaranteed by the definition of the relative sea level, which expresses the water height relative to the ground surface and the variations relative to present-day sea level, \(h_\mathrm{RSL}({ \Omega }, t_\mathrm{pt})\) (1).

One advantage of the spectral representation for the solution is that we are not restricted to the Gauss–Legendre grid on which the sea-level equation is solved, but we can predict the palaeo topography on an arbitrary grid by applying Eqs. 1 and 2. Here, we determine the palaeo topography with respect to the reference topography, \(h_\mathrm{topo}(t_\mathrm{pt})\), according to the resolution of ETOPO1 [2].

The equivalent sea level, which approximates the eustatic sea level according to Mitrovica and Milne [20], is defined as the amount of excess water which is stored during the glaciation in the ice sheets as mass \(m_\mathrm{ice}\) and is expressed as the water column at the present-day ocean surface \(A_O\) with density of water \(\rho _\mathrm{oce}=1020\;\mathrm{kg/m}^3\):

$$ h_{\mathrm{eq}} (t) = -\frac{m_{\mathrm{ice}} (t)} {\rho_{\mathrm{oce}} A(t_{\mathrm{pt}})} $$
(3)

This water level represents a first estimation of the sea-level change during the glacial cycle. We use \(h_\mathrm{eq}\) as eustatic sea level in the following.

To analyse the glacio- and hydro-isostasy predicted by the considered ensemble of earth models, we use the point-wise mean value of the model predictions, \((h^i_\mathrm{RSL} (\Omega, t),\, i=1,\ldots ,70)\). The superscript denotes here the ith ensemble member, the standard deviation, and the range defined as the difference between minimum and maximum prediction. These quantities directly map into \(h_\mathrm{topo}\) (2).

Results and discussion

Fig. 2
figure 2

Variability of relative sea level, \(h_\mathrm{RSL}\), on Northern Hemisphere at last glacial maximum. Left panel the ensemble mean; right panel the corresponding standard deviations. Values below and above the ranges indicated by the given scales are saturated. In both panels, white line is the isoline of mean \(h_\mathrm{RSL}=-100\;m\), dark grey line is the present shore line and orange box marks the studied area

Full size image

Figure 2 shows the distribution of the relative sea level for the Northern Hemisphere at the LGM (20 ka BP, before present). To discuss its spatial pattern, we extended this quantity also over land areas. In regions of positive relative sea level, which coincide with the glacial load distribution (Fig. 1), the GIA process due to ice loading dominates. Here, the lithosphere is displaced downwards due to the acting glacial load (red areas in the left panel). The adjustment is controlled by the viscoelastic earth structure visible in the large standard deviations which appear here due to the considered earth model ensemble used in this study. Away from these regions, the relative sea level is near to −120 m in the ocean basins which agrees with the eustatic sea level proposed by Fleming et al. [10]. When crossing the continental coastlines, the relative sea level rises to values above −100 m as a consequence of hydro-isostasy. Furthermore, we observe an increased standard deviation on the land side.

Considering the −100 m contour in Fig. 2 within the orange box, which marks the region of the Laptev and the East Siberian seas, this region might be dominated by hydro-isostasy. From the right panel, it is visible that GIA has still some influence, as the larger standard deviations at and around the formerly glaciated regions spread into the marked region. The enhanced deviation at the Arctic Ocean is in contrast with other continental coasts where the ocean side shows a smaller standard deviation than the land side.

Fig. 3
figure 3

Variability of palaeo shoreline in studied area at epochs shown at the top left: Green/blue distinguishes land/water, reddish area in between shows variability of predicted shoreline due to earth model ensemble. The color corresponds to that of the isolines representing the standard deviation of relative sea level for the earth model ensemble. Light grey line the present shoreline; dark red line the paleo-shoreline of the ensemble mean and blue line the shore-line prediction based on the sea-level equivalent

Full size image

Figure 3 represents the mean shoreline location and its variability at 20, 15, 11.5 and 5 ka BP, respectively. At 20 ka BP, most parts of the continental shelf were land. The transgression of the shelf seas started around 14 ka BP and was completed around 5 ka BP. This is in accordance with Bauch et al. [5]. The variability in the predicted shoreline location is mainly evident during the transgression process. The third panel shows the distribution at 11.5 ka BP, where the reddish areas show the spread in shoreline location due to the considered earth model ensemble. The area is significantly larger in the East Siberian Sea, whereas the Laptev Sea shows only a small spread. The shoreline based on the eustatic sea level is shown in blue, and does not markedly differ from the predicted shorelines considering GIA—a fact which we will discuss further down.

Fig. 4
figure 4

Relative sea level vs. bathymetry along Profiles A, B and C shown in the top panel. Ensemble mean of relative sea level along the profiles and for denoted epochs is shown in red together with the ensemble range. Present-day bathymetry is shown in grey. The sea-level equivalent is shown as horizontal blue lines for the respective epochs repeated on the right

Full size image

The reason for the different behaviour of the Laptev Sea and the East Siberian Sea can be explained by their different bathymetry. This is shown in Fig. 4 at three S–N profiles crossing the continental shelf in the Laptev Sea (Profile A) at the New Siberian Islands (Profile B) and in the East Siberian Sea (Profile C) where the grey line denotes the respective present-day bathymetry. The profiles were positioned in a way that the present shoreline is located roughly at 200 km distance on each profile. In addition, the figure shows the modelled relative sea level at the four considered epochs and their respective range in red. Also shown is the sea-level equivalent, \(h_\mathrm{eq}\), in blue for the four considered epochs. As \(h_\mathrm{eq}\) is monotonously rising from 20 ka to present, the order follows from top to bottom that for \(h_\mathrm{RSL}\). The locations in the diagram, where the relative sea-level curves cross the present-day topography, denote the coastline at the respective epoch (see (2)).

Profile A shows a rather steep slope on the continental shelf in the Laptev Sea in comparison to Profile C in the East Siberian Sea. In consequence, the ensemble variability does not introduce a wide spread of the predicted shoreline along the profile. The transgression starts here after 15 ka BP. In contrast, the rather shallow East Siberian Sea expanding over 700 km on this profile results in a much wider spread of the predicted relative sea level which is especially evident at 11.5 ka BP. This is emphasised in Fig. 5, where we directly compare the spread of the predicted shoreline position from 20 ka BP to present day for profiles A and C. From this figure, we can quantify the extremely wide spread of predicted shoreline positions to several hundreds of km.

Fig. 5
figure 5

Evolution of shoreline position along Profile A (red) and C (blue). Each line pair represents the shoreline position for the minimum and maximum predicted shoreline of the respective ensemble. Grey line mainly located between the pair of lines denotes the respective shoreline determined from the eustatic sea level

Full size image

The predicted relative sea level shown in Fig. 4 shows a systematic rise towards the coast, which amounts to 15–20 m at 20 ka and remains significant during the transgression phase. This behaviour can be explained by the bending of the lithospheric plate due to hydro-isostasy. Furthermore, the bending starts more offshore in case of the East Siberian Sea (Profile C).

This spatial distribution in predicted sea level is shown in Fig. 6, where the ensemble mean of relative sea level is shown for the same epochs as in Figs 3 and 4. Here, at 20 ka BP, the −100 m isoline is located near the palaeo shoreline at that time and, so, corresponds to the pattern already discussed in Fig. 2. During the transgression phase, this general pattern remains with smaller amplitude.

Fig. 6
figure 6

Ensemble mean of relative sea level in studied area at the considered four epochs

Full size image

The influence of GIA is visible as the higher variability of relative sea level in the offshore part of the profiles of Fig. 4, which cannot be explained by hydro-isostasy and coincides with the adjacent glacial ice sheets of the Kara Sea and Laurentide, respectively (see Fig. 2).

A further interesting aspect shown in Fig. 5 is the coincidence of the sea-level equivalent, \(h_\mathrm{eq}\), in grey, with the location of the predicted shoreline. Although some variability due to the earth model ensemble exists, this fact is surprising, as it shows that the assumption of Bauch et al. [5] to use the eustatic sea level to estimate the transgression history in the Laptev Sea is reasonable. To generate sea-level diagrams, we chose the three coring sites discussed in Bauch et al. [4], and one site at the settlement of Tiksi near the present coast line (Fig. 7). All curves clearly follow the pattern of the sea-level equivalent according to the considered termination described by ICE-5G. With distance from the shelf, we observe the systematic deviation from the sea-level equivalent which amounts to 15 m at the shelf margin (PS2725) and 30 m at the present coast line (Tiksi). Such a deviation of the predicted water depth away from the apparent shoreline is also visible in Fig. 4. This systematic deviation reduces over time, and the sea-level equivalent reaches the \(h_\mathrm{RSL}\) curves at about 14 ka BP. The respective predicted range in relative sea level of the ensemble is shown in the lower panel of Fig. 7. For better visibility, we shifted the curves by 20 m against each other. From this panel, it is evident that the predicted range remains significant until 4 ka BP. Furthermore, it confirms the influence of GIA at the ocean sites, here more visible in the increased variability with distance from the coast, which is due to the influence of the Laurentide and the Kara Sea ice sheet; see also the higher variability in Fig. 3 to the north east and north west.

Fig. 7
figure 7

Relative sea-level curves from 20 ka PB to present at four sites which are shown in the top map. Middle panel mean relative sea-level curves for the sites with colors from the top panel. Thick grey curve the sea-level equivalent due to (3). Lower panel the same curves as in the middle panel, shifted by 20 m from each other to visualise the respective range of the ensemble

Full size image

Conclusions

In this study, we discussed the influence of GIA on the transgressional evolution of the Laptev and the East Siberian seas. Former reconstructions of the transgressional evolution of the these seas were based on the eustatic sea-level change to reconstruct the evolution of the relative sea-level history, whereas it was assumed that this region is located far enough from Scandinavia and North America to be influenced by the ice loading processes of the respective glacial ice sheets. Solving for geodynamically consistent GIA with a spherical symmetric earth structure, where the sea-level redistribution due to mass conservation and change of the geoid is considered, we tested this assumption and reconstructed the sea-level evolution and transgressional history of the East Siberian Shelf region. To allow for uncertainties in the earth structure, we considered a range of earth models by systematically changing the lithosphere thickness as the upper and lower mantle viscosities.

We distinguished two mechanisms in GIA which alter the relative sea level from its eustatic value during time. The first mechanism, glacio-isostasy, is the direct loading response to the changes of the glacial ice sheet. This signal strongly reduces with distance from the glaciated regions and shows a strong sensitivity to the considered earth structure. The second process, hydro-isostasy, is the loading response to the lateral water-load contrast around the continents and results in a flexural response of the lithosphere at the continental margins. In consequence, the sea level towards the land area is altered systematically and the eustatic sea-level drop of 120 m during the LGM results in a sea-level drop of less than 100 m on land.

Analysing these two mechanisms when predicting the sea-level history of the East Siberian shelf, we found that the transgression is definitely dominated by the eustatic sea level with deviations due to hydro-isostasy appearing in sea-level height predictions away from the shoreline position. The surrounding ice sheets influence the signal mainly in the Eurasian basin. Due to the shallow waters on the shelf, the transgression starts around 14 ka BP directly after melt-water pulse 1A [3] and ends around 5 ka BP.

The different continental slopes, a more steep one in the Laptev Sea and a rather shallow one in the East Siberian Sea result in a different response of these shelf regions to the sea-level rise. The flexure of the lithosphere starts more offshore in case of the East Siberian Sea, and, considering different earth structures, we get a wide spread of predicted shorelines during the transgression phase. In the Laptev Sea, the continental slope is steep enough to constrain the shoreline evolution to follow the ensemble mean of the considered earth models with only small spatial variability.

The results show that the assumption to use the eustatic sea level is reasonable to reconstruct the shoreline evolution in the Laptev Sea, whereas the shoreline migration in the East Siberian Sea is strongly influenced by the earth structure considered. Towards the coast, the predicted relative sea level is, due to hydro-isostasy, systematically higher than the eustatic value, and systematically lower towards the ocean, which amounts to 10–15 m difference. At times before 15 ka BP, this coincidence of the predicted paleo-shoreline fails, as the eustatic sea-level estimate is significantly lower than the relative sea level predicted by GIA modelling in the whole region. But due to the steepness of the continental shelf, this deviation has no influence on the shoreline evolution.

The results suggest to use GIA models to analyse the evolution of submarine permafrost, which depends on local variations of sea level and shoreline migrations, in continental shelf seas like the Laptev and East Siberian seas. This dependency is especially of importance when extending the study to the early phases of the Weichselian glaciation, where the Kara Sea ice sheet was much larger [21], and so deviations from eustatic sea level in the Laptev Sea were substantial, and in former glaciation periods when an East Siberian ice sheet is proposed [15, 23].