The effect of sediment loading in Fennoscandia and the Barents Sea during the last glacial cycle on glacial isostatic adjustment observations

. Models for glacial isostatic adjustment (GIA) routinely include the effects of meltwater redistribution and changes in topography and coastlines. Since the sediment transport related to the dynamics of ice sheets may be comparable to that of sea level rise in terms of surface pressure, the loading effect of sediment deposition could cause mea-surable ongoing viscous readjustment. Here, we study the loading effect of glacially induced sediment redistribution (GISR) related to the Weichselian ice sheet in Fennoscandia and the Barents Sea. The surface loading effect and its effect on the gravitational potential is modeled by including changes in sediment thickness in the sea level equation following the method of Dalca et al. (2013). Sediment displacement estimates are estimated in two different ways: (i) from a compilation of studies on local features (trough mouth fans, large-scale failures, and basin ﬂux) and (ii) from output of a coupled ice–sediment model. To account for uncertainty in Earth’s rheology, three viscosity proﬁles are used. It is found that sediment transport can lead to changes in relative sea level of up to 2 m in the last 6000 years and larger effects occurring earlier in the deglaciation. This magnitude is below the error level


Introduction
Erosion in glaciated areas can be larger than in non-glaciated regions (Hallet et al., 1996;Amantov et al., 2011 and references therein), and estimates for sediment deposition in glaciated regions vary from millimeters per year to centimeters per year close to glaciers (Elverhøi, 1984;Finlayson, 2012), which is comparable to global changes in relative sea level during the last glacial cycle (Fairbanks, 1989).Similarly to sea level change, sedimentation rates are enhanced during deglaciation when runoff is larger (e.g., Tucker and Slingerland, 1997;Ivins et al., 2007).These changes in surface loading can lead to changes in sea level and the Earth's solid surface during thousands of years because of viscoelastic relaxation of the mantle.This raises the question of whether erosion and sedimentation that is enhanced during deglaciation affects present-day glacial isostatic adjustment (GIA) measurements.The loading effects of meltwater redistribution are routinely included in models of GIA, but the loading effect of sediment transport is not.Of course, total sea level change is a global effect while sediment transport is a more local effect, and total meltwater volume is much larger than the displaced sediment, which would argue for a smaller effect of sediment loading.On the other hand, sediment density is higher than water density, and the effects of Published by Copernicus Publications on behalf of the European Geosciences Union.
W. van der Wal and T. IJpelaar: Effect of sediment loading in Fennoscandia and the Barents Sea sediment transport during the last glacial cycle could influence present-day GIA measurements locally.
Several studies have investigated the viscous response due to variation in past sedimentation rates.Ivins et al. (2007) force their surface loading model with an estimate of postglacial sedimentation rates of 10 mm yr −1 , compared to a background sedimentation rate over a glacial cycle of 1 mm yr −1 .Their modeling predicted present-day subsidence of 1-8 mm yr −1 although a more recent estimate reduces that amount to 0.5 mm yr −1 (Wolstencroft et al., 2014).Viscoelastic relaxation due to sediment deposition in the Indus River basin and Arabian Sea has been shown to cause changes in relative sea level of up to 2 m over 4000 years (Ferrier et al., 2015).The effect is larger when the entire glacial cycle is considered, which is relevant when sea level data near the deltas are used to constrain global meltwater volume (Ferrier et al., 2015).
While the aforementioned studies focused on sediment loading near river deltas far away from glaciated areas, the large amount of sediment transport involved in glacier growth and melt could also induce paleo-sea-level changes and present-day vertical motion near previously glaciated areas.We refer to material displaced by glacier growth and melt as glacially induced sediment redistribution (GISR).When present-day observations are used to infer viscosity or ice thickness in GIA models, those inferences could be biased when GISR is not taken into account.The objective of this paper is to find out what is the effect of GISR during the Weichselian on GIA observables in Fennoscandia including the Barents Sea.The interest in this region stems from the fact that glacigenic sediment transport is large there (Riis and Fjeldskaar, 1992;Dowdeswell et al., 1996), with the last glaciation depositing sediment layers of up to hundreds of meters' thickness (Elverhøi, 1984).Moreover, several observations of sediment deposition are available from which the loading can be quantified (e.g., Dowdeswell, 1996;Taylor et al., 2002).Here, the focus is on present-day uplift and gravity rate of change and paleo-sea-level data, which are routinely used to constrain GIA models.
Models exist which compute the sediment displacement as a result of the movement of glaciers (e.g., Boulton, 1996;de Winter et al., 2012), but since the ice sheet thickness, as in most GIA models, is not a dynamic model component, erosion is not coupled to the changes in ice thickness in this study.Instead, the amount of GISR is derived from the literature on observed sediment deposits and reported output from a coupled ice-sediment model.Sediment being deposited in the ocean will not only induce vertical motion but also displace water and affect the gravity field.To model this effect, we use the methodology of Dalca et al. (2013) to include sediment redistribution in the sea level equation in a selfconsistent way.Dalca et al. (2013) show that ignoring the time-varying ocean load resulting from sediment redistribution can result in errors in relative sea level (RSL) of up to 40 %.The method will be discussed briefly in Sect. 2. Af-ter that, it is explained how different estimates of GISR are created.Next, sea level change, deformation rates, and gravity rates are calculated for the different sediment transport scenarios and conclusions are drawn about the relevance of GISR for interpreting GIA observations.

Method
The loading effect of ice and meltwater is routinely included in GIA models.The so-called sea level equation is solved, which computes the sea level distribution that accompanies a change in ice volume and corresponding changes in the Earth's shape and gravitational potential field (Farrell and Clark, 1976;Mitrovica and Peltier, 1991).The effect of sediments can also be included in the sea level equation, as shown by Dalca et al. (2013).Here, we follow Dalca et al. (2013), Kendall et al. (2005) and references therein.Only the key elements will be repeated here, and some small differences will be pointed out.
Defining the sea level as the difference between the equipotential surface corresponding to sea level and the solid surface, the sea level (SL) is given by where G is the height of the equipotential surface coinciding with the sea level, R is the height of the Earth's crust, H is the thickness of sediments, and I is the thickness of ice masses supported by land.The aim is to compute the changes in sea level as a result of a changes in total surface mass load L, which is defined as the sum of the changes in mass of water, ice, and sediment: where ρ w , ρ I , and ρ H are the respective densities and S is the ocean thickness.Computing the change in sea level SL requires the change in equipotential surface and the solid Earth displacement, which themselves depend on the change in sea level.The solution requires solving an integral equation which is usually done with an iterative approach.To solve the sea level Eq. ( 2), loading changes are discretized at time steps of typically 1000 years.Two aspects need to be included to ensure accurate representation of surface loads.First, a check is performed at each time step j to see whether ice is grounded or not by requiring that the ice starts to float when the pressure exerted by the ice (prescribed by the ice model) is equal to the pressure of the current sea level.Thus, floating occurs when sea level is positive in the absence of ice and Second, ocean-continent margins change with time to account for ice sheets replacing sea and vice versa, as well as accounting for the change in coastline as sea level rises next to a sloped coastline (see Kendall et al., 2005, and references therein).The change in ocean-continent margins depends on the topography, which depends on the sea level change.This requires an iteration over the complete glacial cycle (the "outer" iteration) on top of the iteration to obtain the sea level change for each time step (the "inner" iteration, denoted with index i).
To start the outer iteration over the glacial cycle, the preglacial topography T 0 is assumed to be equal to the present-day topography T p .With this topography, sea level at each time step is computed.To start the inner iteration for each time step, an initial guess for the change in ocean height is given by where h j is the uniform change in ocean height given by mass conservation with the current ocean basin and C j is the ocean function at time t j .δ denotes a change in one time step which is different from the total change denoted by .
Note that the change in sediment thickness is subtracted here because it is included in the definition of sea level.
After computing sea level increments at all time steps, the topography estimate can be improved using the total sea level rise: With the improved preglacial topography, the computation of sea levels can be repeated (the outer iteration) until the preglacial topography reaches convergence.Erosion will also change the topography, but this effect is not included; the effect of erosion on the location of coastlines is smaller than the loading changes in erosion and sediment deposition themselves, which are the main focus of this study.
To compute the change in equipotential surface and solid surface displacement, the Earth's mechanical properties need to be known.Here, we assume the Earth is radially symmetric, incompressible, and deforming according to a Maxwell rheology.For such an earth model, response functions for an impulse load can be computed in the spherical harmonic domain (Peltier, 1974).An efficient solution method presented by Mitrovica and Peltier (1991) solves the sea level (Eq.2) in the spatial domain while computing the response of the solid Earth in the spectral domain.This method requires transformations from the spatial to the spectral domain in which some accuracy is lost.
The effect of sediment redistribution is implemented in the numerical codes for the sea level equation developed by Schotman (2007).A partial benchmark against other numerical solutions of the sea level equations was carried out in Spada et al. (2012).Rotational feedback is also included in the sea level equation following Wu and Peltier (1984) and Milne et al. (1998).The response of the Earth to surface loading for a radially symmetric Earth is computed with the multilayer matrix propagation normal mode method (Vermeersen and Sabadini, 1997) which is benchmarked in Spada et al. (2011).

Model inputs
The computation requires several inputs, such as elastic parameters, the viscosity of the Earth, and the ice and sediment distributions, which are discussed in the following subsections.For the present-day topography, the ETOPO5 dataset is used.The maximum spherical harmonic degree is 256, and the size of the grid of quantities that are provided in the spatial domain, such as topography and surface load, is 256 × 512.

Model inputs: viscosity and ice loading
In this study we consider a laterally homogeneous Earth model and vary the radial viscosity profile.As a reference profile we use VM5a (Peltier and Drummond, 2008), which is an iteration of the VM2 profile that is used in the creation of ICE-5G (Peltier, 2004).As alternative profiles we select profiles that have been shown by Root et al. (2015b) to provide a good fit to sea level data, GPS, and GRACE (Gravity Recovery and Climate Experiment) data in Fennoscandia.That study found two viscosity profiles: one with higher viscosities in the upper and lower mantle and one with lower viscosities.The fact that sediment loading is not taken into account to obtain viscosity profiles in Root et al. (2015b) will have a minor effect on our results given that three very different viscosity profiles are selected to account for uncertainty in viscosity.Out of those sets we select M8-128-150 and M4-16-80, where the first number denotes the upper mantle viscosity in 10 20 Pa s, the second number denotes the lower mantle viscosity in 10 20 Pa s, and the third number denotes the lithosphere thickness in kilometers.The three viscosity profiles are shown in Fig. 1.Note that viscosity is 10 22 Pa s in VM5a just below the lithosphere, from 60 to 100 km depth.
Since we are only interested in the effect of GISR, the exact ice loading history is of less importance and only influences the effect of GISR through the distribution of meltwater possibly replaced by sediment, which is a smaller effect than the sediment loading itself.For ice loading history, the ICE-5G modelv1.2(Peltier, 2004) is selected, which is provided with time steps of 2000 years from 120 kyr before present (BP) to 32 kyr BP, 1000 years from 32 to 17 kyr BP, and 500 years from 17 kyr BP to the present.

Model inputs: sediment distribution
In order to model the loading effect of GISR, it is necessary to know how much sediment is transported, where it has come from, and where it is deposited.Erosion and deposit W. van der Wal and T. IJpelaar: Effect of sediment loading in Fennoscandia and the Barents Sea estimates for all of Scandinavia during the last glacial cycle are not readily available.Therefore, we created a map of sediment deposition from estimates in the literature of sediment volumes transported in local features: (i) trough mouth fans (TMFs), (ii) large-scale failures, and (iii) basin flux.Each of the features will be briefly discussed in the following.
(i) TMFs are places where rapid flowing ice streams at the end of the continental shelf converge and where sediment is deposited off the shelf; see Fig. 2 for the locations.Local observations come from sonar and seismic profiling (Dowdeswell et al., 1996) and coring (Saettem et al., 1992;Laberg and Vorren, 1996;Taylor et al., 2002;Laberg et al., 2012).Other estimates come from modeling based on bathymetry, elevation, and environmental conditions (Siegert and Dowdeswell, 2002), but these are highly dependent on the amount of ice that is believed to have existed in the Barents and Kara seas.Deposition also takes place on the shelf but is likely smaller than deposition off the shelf (Zieba et al., 2016).
The estimates are compiled in Table A1 in Appendix A.
(ii) Large-scale failures represent the sediment that is displaced after the collapse of the slope.The largest of such events related to the Eurasian ice sheet is the Storegga slide.Haflidason et al. (2004Haflidason et al. ( , 2005) ) estimate the Holocene event to have a volume of 2400-3200 km 3 based on sonar scans and sedimentary cores.A compilation of the studies for this and other slides is provided in Table A2 in Appendix A.
(iii) The TMFs and the large-scale failure are examples of local features.However, constant deposition of sediment from the source area to the basins (Fig. B1) also takes place.During glaciations, sediment activity is increased, with thickness changes estimated to be 5 and 2 cm kyr −1 for the periods of 30-10 kyr BP and 10-0 kyr BP, respectively (Taylor et al., 2002), which amounts to a total volume of 97 and 58 km 3 for the Norwegian basin and 485 and 290 km 3 for the Lofoten basin.There are also channels and canyon systems that are not captured by either the large-scale individual features which, nevertheless, provide larger sedimentation rates than the overall basin estimate.Estimates for the Lofoten channel system amount to 35 km 3 (Taylor et al., 2000), which is small enough that they can be neglected compared to the other events.The basin fluxes are given in Table A3 in Appendix A.
Conflicting estimates are stated for GISR volumes in Tables A1-A3 in Appendix A, for example the Bjørnøya trough mouth fan in Table A1.The timing is also uncertain for most events.Therefore, different sets of GISR loads were created, consisting of minimum, maximum, and moderate estimates from the tables, which are labeled Sed1, Sed2, and Sed3, respectively.The time step is chosen to be the average of the time span given, rounded off to the nearest time step of the ice model.A value of 2300 kg m −3 is taken for sediment density, in agreement with Amantov et al. (2011) and measurements for shallow sediments (Zaborska et al., 2008).The uncertainty in the density is likely smaller than the uncertainty in thickness and timing, which are represented here by the different sediment models.

Results
Figure 3 shows the RSL between the Last Glacial Maximum (LGM) and the present for viscosity profile M4-16-80 together with the location of a selection of RSL sites from the Tushingham and Peltier (1992) database.The RSL mainly represents the solid Earth displacement.It can be seen that the largest subsidence caused by glacial sediment transport is offshore, and the largest uplift is between south Norway and Sweden.Hence the effect on coastal RSL data is limited.Another reason why the effect on RSL data is limited is that the GISR causes a difference in RSL that increases over time.Because the sea level of current sites is set to zero, the largest difference occurs earlier in the deglaciation, which coincides with larger error bars in the data, if records are available at all.This can be seen in Fig. 4, which shows the effect on RSL with and without GISR for the ICE-5G model in combination with the M4-16-80 Earth model.Differences are at the level of several meters, which is below or near the error bar.Values for the other sediment and earth models are given in Table 1.It can be seen that the Amantov model results in values that are a factor 2 or 3 larger.
While the relative sea level over time also includes the geoid effect due to removed mass, the present-day uplift rate only represents the viscous readjustment due to past changes in surface loadings.The pattern of uplift rates is shown in Fig. 5 for the Amantov and Sed1 GISR models.In Fig. 5a Table 1.Maximum effect of sediment loading for different GISR estimates.In each cell the three numbers separated by "|" correspond to earth models M4-16-80 | M8-128-150 | VM5a, respectively.The maximum effect on relative sea level measurements is calculated as the maximum effect at any of the six sites of Fig. 3 at the time at which there are measurements (shown with vertical bars in Fig. 4).The uplift rate is interpolated at the GPS sites of the BIFROST network presented in Lidberg et al. (2010), and the maximum is shown.The maximum positive rate of change in gravitational acceleration in Scandinavia is determined in the land area contained in the box with longitudes from 5 to 37 • and latitudes from 55 to 71 • N; see Fig. 6.The maximum positive gravity rate in the Barents Sea is determined in the box with longitudes between 10 and 100 • and latitudes between 71 and 81  surements can be made.To see the effect on observed uplift rates, the second data column in Table 1 shows the maximum effect of the sediment loading at any of the BIFROST sites of Lidberg et al. (2010).GISR is seen to always increase uplift rates because most of the GPS measurement stations are in previously glaciated areas where erosion took place.Thus, when GPS is used to draw conclusions on or validate glacia-tion history (e.g., Kierulf et al., 2014), the contribution of the ice may be overestimated.Finally, the influence on gravity rates is also investigated (Fig. 6), as gravity rates derived from the GRACE satellite mission constrain GIA in Scandinavia (Steffen et al., 2008;van der Wal et al., 2011) and the Barents Sea (Root et al., 2015a).For a comparison with GRACE data, a maximum spherical harmonic degree of 60 is used in the GIA model, which is the same truncation used in many GRACE studies.Similar to the uplift rate signal, Sed1 reflects the signal of the landslides offshore of southern Norway, and the Amantov model has negative gravity rates west of the Barents Sea where sediment is deposited.To evaluate the magnitude of the effect, Table 1 provides the maximum gravity rate that occurs in the areas that are used for GIA studies: Scandinavia and the Barents Sea; see Fig. 6.The values can be judged by comparing them to the GRACE measurements error.The measurement error of the gravity rates derived from GRACE is computed using the method of Wahr et al. (2004), assuming that residuals obtained after fitting a trend and the secular and annual signal to the monthly gravity fields reflect noise.This method was shown to result in a similar error magnitude to calibrated SDs or a full variance-covariance matrix (van der Wal et al., 2010).The measurement error propagated to the trend has a value of 0.016 µGal yr −1 for a 10year GRACE time series from January 2003 to July 2013.In Scandinavia sediment loading results in gravity rates at the level of the measurement error for the Sed/2/3 models for the three Earth models (between 0.002 and 0.014 µGal yr −1 ) and larger effects for the Amantov model (between 0.007 and 0.027 µGal yr −1 ).However, gravity rates derived from GRACE are not only affected by measurement error.Using data from different processing centers and using different correction models results in a larger spread in the gravity rate than the measurement error.Based on a comparison with GPS data, the rms error was estimated to be 1 mm yr −1 in terms of uplift rate (van der Wal et al., 2011), which is around 10 % of the maximum uplift rates.Assuming most of this re- .Gravity rate due to GISR according to the Amantov model for viscosity profile M4-16-80.Boxes around Scandinavia and in the Barents Sea denote the areas for which the maximum gravity rate is given in Table 1.
flects GRACE errors, this means that the gravity rate error is roughly 0.1 µGal yr −1 .In this light, the values in Table 1 appear to be insignificant, but the effect of sediment loading is systematic.Note that the direct attraction of the current rate of sedimentation in the Barents Sea is not included in our computations.The effect of current sediment transport could become relevant for GRACE studies and should be investigated in future studies.
www.solid-earth.net/8/955/2017/Solid Earth, 8, 955-968, 2017 We investigated the effect of sediment transport during the past glaciation in Scandinavia on current GIA observables.The effect on present-day GIA observables is small compared to the effect of ice loads that are displaced during glaciation.Sediment uptake takes place over a large area, and deposition takes place in limited areas mostly confined to the ocean which can lead to a locally higher signal near areas of sediment deposition.It was found that RSL data are not significantly affected by GISR because those data are located near the shore, in between zones of erosion and deposition, and because a large part of the deposition takes place early in the deglaciation when the errors in RSL data are relatively large.For the LGM, the effect is 18 m; at 6 ka, the effect is below 2 m, comparable to what was found for the Indus River basin (Ferrier et al., 2015).
The effect on present-day uplift rate and gravity rates is also limited.Depending on the estimate for sediment transport that is used, the magnitude of GISR loading effects is near the measurement limit: several tenths of millimeters per year uplift rate and several tenths of microgal per year gravity rate.The magnitude is comparable to a recent estimate of subsidence in the Mississippi delta (Yu et al., 2012;Wolstencroft et al., 2014) but somewhat smaller because of the larger sediment deposition area for the Fennoscandian ice sheet and the earlier demise of the ice sheet in the Barents Sea.The magnitude is smaller than a possible reference frame bias in GPS-derived uplift rates (Lidberg et al., 2010), and in the presence of other GIA model uncertainties, several tenths of millimeters a year do not appear to be significant.Nevertheless, the effect is systematic; correcting for the effect of GISR reduces uplift rates and gravity rates in the land area of Fennoscandia and increases gravity rates off-shore of Fennoscandia.Thus, if uplift or gravity rates are used to infer viscosity profiles or ice thickness, those estimates could be biased.Also, a few tenths of millimeters a year are not negligible compared to global average sea level change.
Lateral variations in earth properties could affect the conclusions.The eastern part of Scandinavia is part of a craton, which manifests itself in large crustal thickness and higher seismic velocities, which extend to the Barents Sea as seen in seismic measurements (see, e.g., Schaeffer and Lebedev, 2013).The large seismic velocities likely result from high viscosity underneath eastern Scandinavia and the Barents Sea (e.g., van der Wal et al., 2013), which could reduce the present-day uplift rate in east Scandinavia and the Barents Sea (Kaufmann and Wu, 1998) but could increase the uplift rate west of Norway.
To correct uplift and gravity rates for the effects of GISR it is necessary for more accurate estimates of sediment transport to be made, including variations in sediment density (Blum et al., 2008;Ferrier et al., 2015), or for ice loading histories currently used in GIA models to be coupled with models for erosion and sediment processes.We suggest investigating the effect of GISR in other areas where last glacial ice caps were located close to the continental shelf and GISR is expected to be large, such as Antarctica and Alaska.

Figure 3 .Figure 4 .
Figure 3. Colors denote the difference between RSL computed for the LGM (26 kyr BP in the ICE-5G model) and computed for the present caused by GISR for sediment models Amantov (a) and Sed1 (b).Note the different color scales.Numbered black dots show the locations of the relative sea level data used in Fig. 7. Earth model M4-16-80 is used for both cases.

Figure 5 .
Figure 5. Uplift rates caused by sediment redistribution according to GISR models Amantov (a) and Sed1 (b).Dots denote locations of the GPS sites of Lidberg et al. (2010).Note the different color scales.Earth model M4-16-80 is used for both cases.

Figure B2 .
Figure B2.Points captured in QGIS from Amantov et al. (2011) with the value of the sediment load expressed in meters of water load.Interpolation is done using splines.

Table A1 .
Estimates for the volume of sediment displaced in different events, with their time span and a reference.No estimates for Bellsund and Kongsfjorden are available, but sediment transport there is expected to be relatively small.

Table A2 .
Estimates for the volume of sediment displaced by largescale failures.

Table A4 .
Sediment estimates created from observations and estimates of individual estimates in Tables A1-A3, as described in the main text.Sed2 and Sed3 are the maximum and minimum models, respectively.