GRACE constraints on Earth rheology of the Barents Sea and Fennoscandia

. The Barents Sea is situated on a continental margin and was home to a large ice sheet at the Last Glacial Maximum. Studying the solid Earth response to the removal of this ice sheet (Glacial Isostatic Adjustment, GIA) can give in-sight to the sub-surface rheology of this region. However, because the region is currently covered by ocean, uplift measurements from the center of the former ice sheet are not available. The Gravity Recovery and Climate Experiment (GRACE) gravity data has been shown to be able to constrain GIA. Here we analyze GRACE data for the period 2003 − 2015 in the Barents Sea and use it to 5 constrain GIA models for the region. We study the effect of uncertainty in non-tidal ocean mass models that are used to correct GRACE data and ﬁnd that it should be taken into account when studying solid Earth signals in oceanic areas from GRACE. We compare GRACE derived gravity disturbance rates with GIA model predictions for different ice deglaciation chronologies of the last glacial cycle and ﬁnd that best ﬁtting models have an upper mantle viscosity equal or higher than 3 · 10 20 Pa · s . Following a similar procedure for Fennoscandia we ﬁnd that the preferred upper mantle viscosity there is a factor 2 larger than 10 in the Barents Sea for a range of lithospheric thickness values. This factor is shown to be consistent with the ratio of viscosities derived for both regions from global seismic models. The viscosity difference can serve as constraint for geodynamic models of the area.

The rheology of the Barents Sea region is expected to be different from that of Fennoscandia, as it borders passive oceanic margins in the north and the west. Seismic tomography reveals lower seismic velocities in Barents Sea than below Fennoscandia (Levshin et al., 2007;Schaeffer and Lebedev, 2013), but not for all seismic periods and depths. 3D viscosity has been implemented in GIA models for the regions and has been found to affect sea level and uplift rates (Kaufmann and Wu, 1998). However, the difference in properties between Fennoscandia and Barents Sea has not been studied explicitly. 5 Constraints from palaeoshoreline data on 1D GIA models resulted in best fitting upper mantle viscosities of 2−6·10 20 P a·s in the Barents Sea region (Steffen and Kaufmann, 2005), while recent work based on RSL data find that best fitting upper mantle viscosity in the Barents Sea region is above 2 · 10 20 P a · s (Auriac et al., 2016). For Fennoscandia, the best fitting upper mantle viscosity is found to be between 3 − 7 · 10 20 P a · s based on RSL data and relaxation time spectra, while best fitting models based on GPS uplift rate measurements have upper mantle viscosities up to 15 · 10 20 P a · s, see the overview in 10 Steffen and Wu (2011). More recent work summarized in Simon et al. (2018) shows an upper mantle viscosity in the range of 3.4 − 20 · 10 20 P a · s. Note that the lower bound for upper mantle viscosity in the Barents Sea is somewhat below that in Fennoscandia. Steffen and Kaufmann (2005) computed RSL misfit and find similar upper mantle viscosity for the Barents Sea and the Scandinavian mainland, but smaller lower mantle viscosity. However, the different studies used different ice histories and relied on multiple data sources, with substantially less coverage in the Barents Sea region. Therefore it is unknown if it can 15 be concluded from previous 1D studies whether viscosity is indeed lower in the Barents Sea than in Fennoscandia.
In this study we analyze GRACE data in the Barents Sea region and Fennoscandia to obtain the GIA signal there, focusing on the first region where the signal to noise ratio is lower. We compare the estimated signal with 1D GIA model output to infer upper or lower bounds in viscosity for different ice deglaciation chronologies. From comparison between the best fitting models for the two regions we draw conclusions on the variation in Earth rheology between the Barents Sea and Fennoscandia.

GRACE Data Processing
Temporal variations of the Earth's gravity field measured by GRACE are related to mass transport within the Earth system due to different geophysical processes, such as hydrology, ongoing cryospheric mass changes, GIA and (post-) seismic signals (e.g., Wouters et al., 2014). To study GIA, other geophysical signals that mask the GIA signal should be removed. Additionally, 25 GRACE data is affected by instrumental noise and the anisotropic sampling of the signal due to the satellites' orbit (Wahr, 2007;Flechtner et al., 2016). Different data-processing techniques have been developed to increase the GRACE signal-to-noise ratio (e.g., Han et al., 2005;Swenson and Wahr, 2006;Kusche et al., 2009). In the following, we detail the post-processing used to analyze GRACE data in the Barents Sea and Fennoscandia with focus on the Barents Sea as it presents additional difficulties due to the smaller magnitude of the signal. 30 In our analysis we use the University of Texas Center for Space Research (UTCSR) release 5 (RL05) (Bettadpur, 2012) up to spherical harmonic degree 60. We use data for the 2003 − 2015 period. We substitute the degree two coefficients with those obtained from satellite laser ranging (Cheng et al., 2013). To increase the signal-to-noise ratio in the Barents Sea we follow the strategy of Root et al. (2015a). We use a Gaussian filter to filter-out the noisy short wave-length gravity data and reduce GRACE's correlated errors which are evident as north-south stripes in the Level 2 data. We also use a high-pass filter in the Barents Sea to remove the long-wavelength signal that contains unmodelled long wave-length phenomena such as global sea-level rise. We adopt the half-widths used in Root et al. (2015a) which were tuned to optimize the signal-to-noise ratio in the Barents Sea (see Root et al. (2015a) Supplementary Material). The low-pass filter half-width ranges from 200 to 300 km while 5 the high-pass filter half-width ranges from 500 to 700 km. As the signal in Fennoscandia is larger and has a larger wavelength we only use the low-pass filter there with a half-width also ranging from 200 to 300 km.
We compute gravity disturbance rate (gravity rate in the following) (Hofmann-Wellenhof and Moritz, 2006). We use the least square method to obtain the secular, annual and semiannual signals of each time series of Stokes' coefficients. We estimate GRACE measurement errors (σ GRACE ) using the residuals after the secular, annual and semiannual signals are removed from 10 the signal .
After processing the signal as explained above, the GIA signal is evident as a positive gravity rate in Fennoscandia and the Barents Sea ( Figure 1). However, this signal cannot be directly interpreted as it contains the trend of other geophysical processes as well, one of them being hydrology. Secular changes in land water storage result in gravity trends that should be subtracted when analyzing GRACE data in continental areas. The long-term hydrology signal in Fennoscandia is probably 15 small, as demonstrated by the good agreement between GIA signal derived from GRACE and GPS (van der Wal et al., 2011).
However, the hydrology signal of the Russian Arctic Archipelago (Novaya Zemlya, Franz Josef Land and Severnaya Zemlya) can leak into oceanic areas. We subtract the hydrology signal using the GLDAS hydrology model (Rodell et al., 2004). Because its reliability for the islands of the Arctic Archipelago is not well-known, we follow Matsuo and Heki (2013) and take the amplitude of its trend in the Barents Sea as an indication of the uncertainty in the hydrology signal in these polar regions 20 (σ hydrology ).
Present-day changes in the cryosphere and the resulting present-day solid Earth response can also mask the GIA signal. In particular, the glaciers of the islands Svalbard and the Russian Arctic Archipelago are experiencing significant mass changes evident in GRACE observations which partly mask the GIA signal in the Barents Sea region (see Figure 1). Independent data on mass changes in Svalbard and the Russian Arctic Archipelago is limited. Moholdt et al. (2012) derived trends using ICESat 25 for the 2003-2009 period using altimetry; other authors (e.g., Schrama et al., 2014;Matsuo and Heki, 2013) have used GRACE data. For the period 2003-2008, GRACE estimates are lower than altimetry estimates but agree within uncertainty (Root et al., 2015a). In Simon et al. (2018) ice mass loss estimates from altimetry and glaciology for a longer period were shown to be much larger than GRACE estimates in Svalbard, Franz Josef Land and Novaya Zemlya, and the former were scaled down in that study. Here we follow Root et al. (2015a) and use ice loss corrections obtained using the mascon method of Schrama et al. 30 (2014) (see Table 1) to remove the ice loss signal taking into account elastic loading (Wahr et al., 1998).
To obtain the present-day mass changes from GRACE, a GIA correction needs to be first applied As our aim is to quantify the GIA signal in the central Barents Sea, the problem seems circular. However, the GIA model has a relatively small effect on the derived present-day mass changesWe account for uncertainity in mass loss estimations due to GIA by employing an ensemble of ice deglaciation chronologies and Earth rheological parameters. We use the ICE-5G model and two runs of the 35 GSM (Tarasov et al., 2012) with maximum and minimum ice sheet extents combined with the VM5a Earth model (Peltier, 2004) and an Earth model with a stronger mantle, as well as the W12 ice model (Whitehouse et al., 2012) with a strong mantle. Mass loss changes obtained using the different GIA models are shown in Table 1: more massive ice sheet models and stronger mantles result in higher mass loss rates. The error in the derived mass changes due to uncertainty in GIA is similar to the GRACE measurement error. We use the error bars of the estimated mass changes for Svalbard and the Russian Arctic 5 Archipelago to estimate the error in the recovered GIA gravity rates due to uncertainty in mass loss changes in the region (σ ice ).
Finally, for the Barents Sea, the Greenland mass loss is already filtered when using the high-pass filter, but for Fennoscandia we need to remove it. To do so we use ICESat mass changes from Sørensen et al. (2011).
We account for the uncertainty in non-tidal ocean changes by using the ECCO ocean model (Forget et al., 2015) as alternative for the ocean model used in standard GRACE level 2 processing. In that case we first add back the GAB products to restore 10 the full GRACE ocean mass signal (Flechtner et al., 2015;Yu et al., 2018) before subtracting the ECCO ocean model. The ECCO model is a dynamically consistent ocean model constrained with observations from altimetry, Argo floats and GRACE.
The model has been shown to correctly capture long-term bottom pressure variability in the Arctic Ocean and Adjacent Seas (Peralta-Ferriz, 2012). The version of the ocean model we use is the ECCOv4-llc270 compilation. This compilation covers the period 2001-2015 which means the GRACE time-series that we use in the Barents Sea is limited to this period. We obtain 15 gravity rates in the central Barents Sea using the UTCSR GRACE solution corrected with both the OMCT and the ECCO ocean models. The differences between these two solutions are used as an indication of the uncertainty in non-tidal ocean changes (σ ocean ).
We estimate the total error in the gravity trends by assuming that the different error sources are uncorrelated: (1) 20 The assumption that errors are uncorrelated requires further discussion. GRACE data is assimilated in the ECCO ocean model.
However, GRACE is only one of the 40 data sets used in the inversion process and the final product does not fit GRACE data well (Yu et al. 2018). Therefore there will be only a weak correlation with the GRACE data used in our estimation. Correlation between land surface hydrology models and present-day ice melt is not expected, because hydrology models have little skill in predicting trends and do not model areas of permanent snow. Finally, ice loss changes errors (σ ice ) arise due to uncertainty in 25 the GIA model and GRACE measurement error, we cannot rule out that the second error component might be correlated with For the Barents Sea we consider the four terms; while for Fennoscandia we only consider GRACE measurement errors, as the ice loss changes in the Arctic Archipelago and ocean bottom pressure changes have a very small effect on the gravity trends recovered in Fennoscandia, and Greenland's mass loss signal is well known from altimetry measurements.

GIA Modelling
We compare GRACE derived gravity rates with those predicted by GIA models. To compute the gravity trends the sea level equation is solved, using the pseudo-spectral approach presented in Mitrovica and Peltier (1991). We use the same code as Barletta and Bordoni (2013). To be able to run calculations for many different Earth parameters and ice models we assume that solid Earth properties only vary radially, which allows to compute GIA response for different regions separately with different 5 viscosity profiles, but neglects effects of viscosity changes in surrounding regions. While this approach has been used in other GIA studies (Lambeck et al., 1998;Steffen et al., 2014), it has been suggested that far-field viscosity variations are relevant in Fennoscandia (Whitehouse et al., 2006).
We neglect the loading effect due to sediment transport during deglaciation, as the effect is small and well below that of the unknown ice thickness (0.01 to 0.05 µGal/yr in Fennoscandia, and below 0.014 µGal/yr in the Barents Sea as shown in 10 van der Wal and IJpelaar (2017) ). To study the effect of the ice deglaciation history on the present gravity rates we start by using a reference Earth model based on the averaged VM2 model which is similar to the VM5a model (Peltier, 2004;Argus et al., 2014). The model consists of a 90 km lithosphere, a 570 km upper mantle with a viscosity of 0.5 · 10 21 Pa · s and a 2216 km lower mantle with an average viscosity of 2.6 · 10 21 Pa · s. The elastic properties of the Earth are based on the PREM model (Dziewonski and Anderson, 1981). To investigate the effect of the Earth's rheology, we vary the upper mantle viscosity 15 between 0.1−1.6·10 21 Pa · s and the lithospheric thickness between 40−180 km (Table 2). We do not change the lower mantle as its viscosity cannot be constrained uniquely from data in Fennoscandia (Steffen and Kaufmann, 2005). However, gravity rates are influenced by the lower mantle viscosity which is discussed in the comparison with other GIA studies for the region.
We use an ensemble of ice histories that reflects the uncertainty in the deglaciation history of the European Ice Sheet Complex (EISC), the amount of ice in the SBKIS and the Scandinavian Ice Sheet (SIS) for the different ice deglaciation scenarios is 20 shown in Figure 2 . The ice sheet models that we use can be divided in two main categories: (1) empirical ice sheet models based on GIA observables and empirically-determined ice extents, and (2) those based on numerical ice-sheet modeling forced under different palaeo-climate scenarios and tuned to fit different constrains. A fundamental difference between these two kinds of models is that GIA-based palaeo-ice sheet models are explicitly associated with a specific Earth model. The first set of models is represented by the ICE-5G and ICE-6G models (Peltier, 2004;Peltier et al., 2015;Argus et al., 2014). Both models 25 start the ice build-up 122 ka BP. The second set consists of three models obtained using the Glacial System Model (GSM) for Northern Europe (Tarasov et al., 2012), the University of Tromsø Ice Sheet Model (UiT ISM) (Patton et al., , 2017, and the S04 ice sheet model (Siegert and Dowdeswell, 2004), which are further described below.
The three ice sheet models obtained using the GSM model are a subset of a bigger ensemble used in Root et al. (2015a) which showed good agreement with GRACE observations. The ensemble was obtained using a Bayesian calibration of GSM 30 runs with RSL curves, present-day ground velocities and ice deglaciation margins from the DATED-1 project (Hughes et al., 2016). The VM5a rheology model was employed as reference during the calibration process, however, errors introduced by the rheology model were accounted for during the calibration process which implies that this model is not as strongly biased by a single viscosity profile as the ICE-5G and ICE-6G models. The three selected models consist of a late deglaciation model, labeled nn45283 in Root et al. (2015a) and two early deglaciation models, nn56536 and nn56597, with different maximum ice volumes. The build-up phase is faster than for the ICE-5G and ICE-6G models, build-up starts 28 ka BP. The models will be labeled T1 (nn45283), T2 (nn56536) and T3 (nn56597) to simplify the notation.
The University of Tromsø Ice Sheet Model is based on a 3D thermomechanical ice model which uses an approximation of the Stokes equations forced by climatic and eustatic sea level perturbations to simulate the evolution of the EISC. The 5 model is constrained using different geophysical and geological data sets including geomorphological flow sets, moraine and grounding zone wedge positions and isostasy patterns and is consistent with the DATED-1 ice sheet margins. Isostatic loading is implemented using the elastic lithosphere/ relaxed astenosphere model of Le Meur and Huybrechts (1996). The model has no ice in the region before 37 ka BP.
Finally, we consider an ice sheet model which gives a lower bound for the mass present in the Barents Sea during the LGM, 10 the S04 model (Siegert and Dowdeswell, 2004). The model is based on the continuity flow equations coupled with a model of water, basal sediment deformation and transportation. The model is forced with eustatic sea level curves of the last 30 ka and palaeo air temperatures and precipitation and assumes an ice-free scenario before 32 ka BP. Bedrock topography is adjusted for isostasy using the method of Oerlemans and van der Veen (1984).
The only global ice sheet models are the ICE-5G and ICE-6G, for the other ice sheet models we use the ICE-6G ice model 15 outside the EISC. We include the build-up and deglaciation phase of the last glacial cycle. All ice sheet models are sampled in a grid with a spatial resolution corresponding to a 128 degree Gaussian grid and the output of the model is truncated at degree 60 and processed using the same filters used to process the GRACE data.

Model Performance Assessment
We assess the fit of the modelled and estimated gravity rates for different combinations of ice deglaciation history and rheology. 20 As GRACE's resolution is of the same order of magnitude as the extension of the SBKIS we cannot resolve the differences in the shape of the ice sheet in the data. Thus we assess the model fit only by comparing the maximum modelled (m i ) and estimated (e i ) gravity rate in the central Barents Sea and Fennoscandia and normalize this difference using the observation error (σ i ). In order to make the results as independent of the filter parameters as possible, we compute the average of the misfit obtained using different filter configurations: where N is the number of filter settings. For each ice sheet model we define a 95% confidence interval for the solid Earth parameters as all models with a ∆χ 2 smaller than 5 (Press et al., 1992): The low pass-filter is varied between 200 and 300 km half-width in 20 km intervals. Additionally, in the Barents Sea the high-pass filter is varied between the 500 and 700 km half-width in 100 km intervals.

GRACE GIA signal in Fennoscandia and the Barents Sea
We use the methods presented in Section 2 to obtain the gravity rates over Fennoscandia and the Barents Sea. A clear positive 5 anomaly is evident both in Fennoscandia and the central Barents Sea where the main domes of the Scandinavian Ice Sheet and Svalbard-Barents-Kara Ice Sheet were presumably located (Figure 1), we assume this signal to be entirely due to GIA and call it the estimated GIA signal. The melting of ice in Svalbard and the islands of the Russian Arctic Archipelago is also evident as a negative gravity trend. After removing the mass loss signal as explained in Section 2, we observe that most of the signal of Novaya Zemlya, Svalbard and Franz Josef Land is indeed removed (Figure 1). However, there is still a negative gravity rate 10 left over Severnaya Zemlya, indicating that our ice loss changes might be underestimated for this island. We do not observe a clear positive signal in the Kara Sea, which indicates that if it was glaciated during the LGM the amount of ice present there was much smaller than that located in the Barents Sea. This fact advocates against the larger ice sheets in Denton and Hughes (1981); Grosswald (1998); Grosswald and Hughes (2002) and further confirms the results of the DATED-1 (Hughes et al., 2016) and QUEEN projects (Svendsen et al., 2004). 15 We obtain the maximum gravity rate in the Barents Sea for different filter configurations using the OMCT and ECCO ocean models. Figure 3 shows the maximum gravity rates for a 500 km high-pass filter and different low-pass filter half-widths. As expected, we observe that the maximum gravity signal reduces with increasing filter half-width and so does the error. The gravity rates recovered using the ECCO ocean model are systematically higher that those obtained with the OMCT model. We also show a breakdown of the error (Figure 4) for different low-pass filter half-widths. We observe that the hydrology signal 20 leaking into the Barents Sea is very small and the error budget is dominated by the uncertainty in present-day ice changes, the GRACE measurement error and the non-tidal ocean signal. Moreover, we observe that while the other error sources decrease with increasing filter half-width the ocean error does not. This implies that it has a wavelength similar to that of the GIA signal we want to resolve.
3.2 Implications for viscosity and ice sheet chronology 25 We perform three experiments. In the first experiment we only study the effect of the ice history on the model misfit. We use the reference Earth model (see Table 2) and compare the fit of the predicted gravity rates for different ice deglaciation models with the GRACE derived gravity rate. In the second experiment, we change the Earth rheological parameters to obtain the subset of ice deglaciation histories and Earth rheological parameters that best fit the GRACE observations. Thirdly, we repeat the second experiment for Fennoscandia and compare the optimal solid Earth parameters for both regions to detect possible variations in 30 rheological parameters. Figure 3 compares the maximum present-day estimated gravity rates in the Barents Sea with those given by the different ice sheet models. It must be noted, that the maximum gravity rates produced by each ice history are not only related to the maximum ice volume attained during LGM, but also its geographical distribution and the onset of the deglaciation process.
As an example, we find that while the T2 model has more ice in the Barents Sea than the T1 model, it results in lower gravity rates. This is because deglaciation starts earlier in the T2 model than in the T1 model when the sensitivity of the present gravity 5 rates to mass changes is higher as shown in Figure 5. Similarly, the highest gravity rates are associated with the UiT ISM even though it has less ice than the ICE-5G and ICE-6G model. This is because the UiT ice sheet model has more ice in the central Barents Sea during the last phase of deglaciation. In fact, the model includes an ice bridge between Svalbard, Franz Josef Land and Novaya Zemlya with ice thickness as large as 2000 m at 14.5 ka BP which does not disappear until 12 ka BP. This is not present in either the ICE-6G or the ICE-5G models.

10
When we compare the modelled and estimated gravity rates we find that, for the reference Earth model, the T1, T2 and T3 ice sheet models are the closest to observations. The S04 ice sheet model performs worse; the model does not have enough ice in the region. This result is in accordance with Auriac et al. (2016) who found poor agreement between the S04 model and RSL curves. The more massive ICE-5G, ICE-6G and UiT models result in gravity rates that are too high. However, the discrepancy between these models and GRACE derived estimations is reduced if we use the ECCO ocean model instead of the OMCT. 15 Furthermore, GRACE data can be reconciled with the UiT ISM if the maximum volume of mass in the model is reduced by around 1 m of equivalent sea level rise or if deglaciation started 1 kyr earlier.
Next, we study the effects of changing the solid Earth rheology in the Barents Sea. Figures 6 and 7 (left column) show the misfit of the different ice sheet models to the estimated maximum gravity rates in the Barents Sea for different rheology models. We see that there is a large subset of Earth rheological parameters for which the modelled gravity rate is within the 20 95% confidence interval of the GRACE estimated gravity rate.
The T1, T2 and T3 ice sheet models present a good fit to the observations for a large subset of Earth models including the reference Earth model (ν = 5 · 10 20 Pa · s, h = 90 km). For the less massive S04 model the confidence interval extends from ν = 5 · 10 20 to ν = 1.6 · 10 21 Pa · s. In contrast, for the more massive ice sheets (ICE-5G, ICE-6G and UiT ISM) the subset of Earth models which present a good fit to the Barents Sea observations is smaller and does not contain the reference 25 Earth model. These models, however, fit the observations either for a less viscous upper mantle or for a thicker lithosphere when upper mantle viscosity is fixed. If a less viscous upper mantle viscosity is used the relaxation time of the solid Earth is decreased and the sensitivity to mass changes that occurred during the LGM decreases (see Figure 5). On the other hand, a thicker lithosphere acts as a low pass filter that smooths the gravity signal, reducing its maximum value.
Our results for the UiT ISM are consistent with those obtained by Patton et al. (2017) who inferred an upper mantle viscosity 30 of 2 · 10 20 Pa · s based on RSL data. The lower bound obtained with the other ice models is similar as models with a low viscosity have little sensitivity to mass changes during the early deglaciation phase, where differences between ice models are more manifest (Figure 2). Overall, . for a lower mantle viscosity of 2.6 · 10 21 Pa · s, we obtain a lower bound for the upper mantle viscosity of 2 · 10 20 Pa · s, which agrees with the range of possible upper mantle viscosity found in Auriac et al. (2016) using RSL curves and GPS uplift measurements. We refrain from drawing conclusions on the preferred lithosphere thickness 35 from the misfit plots because the lithosphere has a large influence on the shape of the gravity rate pattern which was not used as constraint here. A higher lower mantle viscosity can result in a lower bound as shown in Steffen et al. (2010); Root et al. (2015b).
We follow the same procedure for Fennoscandia to obtain the subset of Earth rheological parameters and ice sheet deglaciation histories with an acceptable agreement with the GRACE observations (see right column of Figures 6 and 7). It must be 5 noted that the values of the χ 2 are higher for Fennoscandia than the Barents Sea and thus the confidence interval is smaller.
The reason is twofold: the observation error is smaller as compared with the Barents Sea, where uncertainty from mass changes in the glaciers of the surrounding islands and non-tidal ocean changes increase the error bars; and the GIA signal is higher in Fennoscandia than in the Barents Sea (see Figure 1). Nevertheless we can compare the best-fitting models for both regions..
We observe that, contrary to what we got for the Barents Sea, the combination of the ice sheet models ICE-5G and UiT with 10 the reference lithospheric thickness and upper mantle viscosity have a good fit fit. (Figures 6 and 7). As already mentioned, the ICE-5G and ICE-6G models have been constrained using GIA observations, which are abundant in Fennoscandia. As we are using these models with an Earth rheology similar to its reference rheology it is not surprising that the ICE-5G model presents a good fitin this region, however the ICE-6G model performs better with a more viscous mantle due to its lower ice volume. The T1-3 models do not fit the estimated GIA signal with the reference Earth model and require a more viscous mantle. The early 15 deglaciation of the S04 model results in low gravity disturbance rates that do not fit the GRACE estimated gravity disturbance rates. For Fennoscandia we find a lower bound for the upper mantle viscosity of 5 · 10 20 Pa · s which is consistent with current estimates (Simon et al., 2018).
We can infer lateral rheology changes by comparing the optimal Earth rheological parameters obtained for both regions. For each ice deglaciation chronology, we compare the two confidence intervals as well as the best fitting upper-mantle viscosity 20 obtained for each lithospheric thickness. We observe that for the UiT, ICE-5G and ICE-6G model both the confidence intervals as well as the best fitting models systematically prefer a less viscous upper upper mantle in the Barents Sea as compared with Fennoscandia. This is also the case for the T1,T2 and T3 models when the best fitting models are compared, although there is an overlap of models of high upper mantle viscosity and thick lithospheres with a good fit in both regions. This systematic difference is likely evidence of lateral variation in Earth rheology.

Lateral viscosity variation
To strengthen the conclusion of viscosity differences between the two regions, we derive viscosity estimates in an independent way, based on seismic velocity anomalies and experimentally derived flow laws. The absolute viscosity values obtained in this way contain large uncertainty, but the relative difference resulting from the seismic models should represent real change in temperature or composition. Therefore we focus on the ratio between the viscosities beneath Fennoscandia and the Barents 30 Sea and check whether it agrees with the outcome of the GIA model misfit.
To take uncertainty in seismic velocity anomalies into account we use two global seismic tomography models: S40RTS Ritsema et al. (2011) and Schaeffer and Lebedev (2013) (labeled SL) which has higher spatial resolution but reduced sensitivity with depth. For both, the reference model is adjusted to account for a jump in the seismic velocity anomaly in PREM (Dziewonski and Anderson, 1981) and AK135 (Kennett et al., 1995) below 200 km. Shear wave velocities are converted to temperature using relations from geochemistry (Goes et al., 2000;Cammarano et al., 2003) for primitive mantle composition and accounting for anelasticity (anelastic correction model Q4 from Cammarano et al. (2003)). Differences in composition between the Barents Sea and Scandinavia could play a role, but is unlikely to reverse the temperature contrast, due to the first order effect of temperature on seismic velocities in the upper mantle (Goes et al., 2000). To compute viscosity we follow the 5 procedure in Wal et al. (2013) and insert temperatures in the olivine flow laws of Hirth and Kohlstedt (2013). The flow laws for diffusion and dislocation are added, which means the viscosity depends on grain-size and stress. Stress is taken from a 3D GIA model which uses the ICE-5G ice load. Background stresses due to mantle convection are neglected as recent work suggest little interaction between GIA and mantle convection (Huang, 2018). Grain size is chosen to be 4 mm or 10 mm. 4 mm gave best overall fit to GIA data in and 10 mm grain size resulted in the best fit with the observed maximum uplift rate (Wal et al., 10 2013).
To be able to compare against viscosity for the upper mantle in the previous section we use viscosity averaged between 225 and 325 km. This depth is a trade-off; shallower layers have lower temperature and small viscous deformation during the glacial cycle, while for deeper layers the seismic models are less accurate. The depth range is also close to the depth to which the gravity rate in Fennoscandia is most sensitive, see the sensitivity kernels in van der Wal et al. (2011). The viscosity maps 15 are plotted in figure 8. In principle all viscosity values around the ice load play a role in the GIA process, but the highest sensitivity is to values directly underneath the ice load (Paulson et al., 2005;Wu, 2006). We compute the average of viscosities for the locations where LGM ice heights are above 1500 m which covers most of the land mass of Scandinavia, and most of the Barents Sea (see dashed brown contour). Viscosity is computed separately for the region below 71 • latitude for Fennoscandia and above for the Barents Sea. We find that the average viscosity below Fennoscandia is a factor of 2.3 to 2.4 times higher than 20 that in the Barents Sea. This agrees well with the change in best fit upper mantle viscosity that can be seen in the misfit figures 6 and 7. There could still be an effect of 3D structure that is not captured by modelling both regions with 1D models, such as lateral variations within Fennoscandia (Steffen et al., 2014) or the influence of viscosity from outside each region.

Conclusions and Discussion
In this study, we analyse GRACE data in the Barents Sea to constrain the Earth rheology in the region. We compare the fit of 25 different GIA models in Fennoscandia with that for the Barents Sea to find if there is a difference in viscosity between the two regions. We investigate several deglaciation chronologies of the SBKIS, some of which are not explicitly tied to a viscosity model. We use GRACE data for the period 2003 − 2015 and process it to reveal the GIA signal. The ice loss signal from the Svalbard and the Russian Arctic Archipelago is removed using mass change values obtained from GRACE using the mascon method. We observe a positive gravity anomaly in the Barents Sea but no significant anomaly in the Kara Sea, which shows 30 that the ice cover at LGM was considerably thinner there than in the Barents Sea, in agreement with recent studies.
The Barents Sea GIA signal is in a region now covered by sea; therefore, the gravity trends might be affected by non-tidal oceanic mass changes. We correct GRACE gravity rates in the Barents Sea using either of two ocean models, the OMCT and ECCO ocean model, and find higher gravity rates using the ECCO model. The difference in the ocean signal according to the two models is large in the Barents Sea. This uncertainty has not been considered in previous studies of the GIA signal in the region (e.g., Root et al., 2015a;Simon et al., 2018;Kachuck and Cathles, 2018) and thus the errors bars in those studies were provably underestimated. This result has also implications for GRACE studies of non-oceanic mass changes, such as post seismic deformations, in ocean areas (e.g., Han and Simons, 2008;Wang et al., 2012) which possibly have higher uncertainty 5 than previously thought due to errors in the ocean model.
We compare the GRACE derived gravity rates with modelled ones to infer geophysical constraints for the Earth rheology and ice sheet chronology in the Barents Sea region. For a three-layer average of the VM2 viscosity profile (Peltier, 2004) we find, as Root et al. (2015a), that thick ice sheet models (ICE-5G, ICE-6G and UiT) do not fit GRACE observations, while the less massive ice models (T1,T2 and T3) do. Upper mantle viscosity and lithospheric thickness was varied for each ice sheet 10 chronology between 0.1 · 10 21 − 32 · 10 21 Pa · s and 40 − 180 km. We find that the ICE-5G, ICE-6G and UiT ice sheet models can be reconciled with GRACE observations provided the upper mantle viscosity is lower or the lithosphere thicker than in the VM2 model. The same conclusion is reached in Auriac et al. (2016) using GPS uplift measurements and RSL curves instead of gravity data.
The interplay between ice deglaciation chronology and Earth rheology makes it difficult to constrain the ice deglaciation 15 chronology in the Barents Sea (Kachuck and Cathles, 2018). Root et al. (2015a) used GRACE data to conclude that the SBKIS had less ice than previously thought (5 − 6.3 m of equivalent sea level versus 8.3 m). To do so, they used ICE-5G and ICE-6G and showed that they do not obtain the estimated gravity rate when these ice models are combined with their corresponding Earth rheology model. However, here we use the UiT ISM which does not come with an a priori Earth rheoloy model and which contains around 7.5 m of equivalent sea level rise and show that it can fit GRACE observations provided the upper 20 mantle viscosity is around 3 10 20 Pa · s if the lithophere is thinner than 130 km. However, we are able to place a constraint on upper mantle viscosity. From the misfit of all investigated ice chronologies and using a lower mantle viscosity of 2.6 · 10 21 Pa · s, we infer a lower bound for the upper mantle viscosity of 2 · 10 20 Pa · s, which agrees with previous constraints derived from RSL and GPS uplift observations Auriac et al. (2016).
We also study the misfit of GRACE observations to the GIA models in Fennoscandia. For a 2.6 · 10 21 Pa · s lower mantle 25 viscosity, we obtain a lower bound of 5 · 10 20 Pa · s for the upper mantle viscosity, which is consistent with current estimates.
Given all the ice sheet deglaciation chronologies we find that the lower bound for the upper mantle viscosity is a factor of two smaller in the Barents Sea (or, alternatively, the lithosphere thickness should be increased there). Unless all the tested ice deglaciation chronologies are biased in the same direction, this result is evidence of lateral changes in viscosity in between the two regions.

30
To strengthen the finding of viscosity difference between the two regions, we compare our results with viscosity derived from global velocity anomalies and flow laws for mantle material and find that the average viscosity in the Barents Sea is a factor of 2.4 lower than in Fennoscandia. This agrees very well with the results derived from the misfit of GIA models to GRACE data, and strengthens the conclusion that there is a small but significant difference in average upper mantle viscosity between the two regions. This findings have implications for ice sheet models inverted with just one viscosity profile (e.g.,

35
ICE-5G, ICE-6G) and advocates in favour of including lateral Earth rheological parameters in GIA models. The constraints on viscosity variations can be also used to calibrate other geodynamic models of the regions.
The GIA signal for different ice deglaciation histories with a reference Earth model is also shown.