Glacial‐Isostatic Adjustment Models Using Geodynamically Constrained 3D Earth Structures

Glacial‐isostatic adjustment (GIA) is the key process controlling relative sea‐level (RSL) and paleo‐topography. The viscoelastic response of the solid Earth is controlled by its viscosity structure. Therefore, the appropriate choice of Earth structure for GIA models is still an important area of research in geodynamics. We construct 18 3D Earth structures that are derived from seismic tomography models and are geodynamically constrained. We consider uncertainties in 3D viscosity structures that arise from variations in the conversion from seismic velocity to temperature variations (factor r) and radial viscosity profiles (RVP). We apply these Earth models to a 3D GIA model, VILMA, to investigate the influence of such structure on RSL predictions. The variabilities in 3D Earth structures and RSL predictions are investigated for globally distributed sites and applied for comparisons with regional 1D models for ice center (North America, Antarctica) and peripheral regions (Central Oregon Coast, San Jorge Gulf). The results from 1D and 3D models reveal substantial influence of lateral viscosity variations on RSL. Depending on time and location, the influence of factor r and/or RVP can be reverse, for example, the same RVP causes lowest RSL in Churchill and largest RSL in Oregon. Regional 1D models representing the structure beneath the ice and 3D models show similar influence of factor r and RVP on RSL prediction. This is not the case for regional 1D models representing the structure beneath peripheral regions indicating the dependence on the 3D Earth structure. The 3D Earth structures of this study are made available.

The main parameters controlling the GIA process are the glaciation history representing the forcing and the Earth structure representing the viscoelastic behavior of the solid Earth. The consideration of geological reconstructions and geodetic observations plays an important role in constraining factors associated with the glacial loading and Earth structure and for the validation of model parameterization (e.g., van der Wal et al., 2015). But substantial uncertainties in the forcing and the parameterization of the response functions as the fading memory of viscoelastic material demand independent constraints on ice history or Earth rheology (Whitehouse, 2018). In particular, seismic tomography and geodynamics provide independent constraints on solid-Earth parameterization.
The majority of GIA models are based on the assumption of solid-Earth structures that vary with depth alone. These one-dimensional (1D) GIA models are applied in regional studies (e.g., James et al., 2000) as well as in global studies. In concert with constraints on the glaciation history, such 1D Earth structures are tuned to derive a GIA model which can fit observational data Kaufmann & Lambeck, 2000;Lambeck et al., 1998Lambeck et al., , 2014Mitrovica & Forte, 1997;. More recent studies focus on ensemble approaches. So, Yousefi et al. (2018Yousefi et al. ( , 2021 focused on the Pacific coast of central North America. They combined 29 ice-sheet reconstructions and more than 700 1D Earth structures to set up a large ensemble, but failed to achieve an acceptable fit to a global set of RSL data with a single set of model parameters, suggesting a requirement for a more complex Earth structure. From seismic tomography models, which contain information about seismic wave speeds, and derived temperature and viscosity structures of the Earth, it is known that the Earth structure is laterally heterogeneous (e.g., Kaban et al., 2007;Schaeffer & Lebedev, 2013;Steinberger & Calderwood, 2006). Variations in viscosity likely exceed three orders of magnitude (Kaban et al., 2007). Compared to 1D structures, 3D structures provide improved representations of the structural features in different regions of the Earth. But, studies on seismic tomography models and 3D parameterization of the Earth interior reveal the existence of several degrees of freedom, for example, for the determination of lithospheric thickness and the mapping from seismic velocity variations to temperature variations and to viscosity variations (Ivins & Sammis, 1995;Trampert & van der Hilst, 2005). A recent study developed an inverse calibration scheme based on experimental results of seismic energy propagation (anelasticity) in polycrystalline materials to map seismic tomography to temperature and viscosity Richards et al., 2020). Furthermore, tomography models do not exactly represent the seismic structure due to limited resolution, and although some features such as fast velocities in continental cratons as well as variations of seismic velocities with ocean floor age consistently appear, there are still considerable differences between recent models (see Hosseini et al., 2018, <https:// www.earth.ox.ac.uk>). A velocity-to-temperature conversion is probably more appropriate in the sub-lithospheric mantle than in the lithosphere, because the latter is likely chemically distinct (Forte & Perry, 2000;Jordan, 1988). Also, chemical heterogeneities are likely reduced in the sub-lithospheric mantle, due to mixing. Lastly, there are also large uncertainties in the temperature dependence of viscosity. In particular, there are several deformation mechanisms, including diffusion and dislocation creep, which contribute to the rheological behavior and result in a nonlinear stress-strain relation, such that viscosity is an effective parameter that depends on the assumed stress state and strain rates.
The importance of lateral variations in Earth structure has been discussed since the 80s (Sabadini et al., 1986), and the impact of a 3D Earth structure on GIA has received significant recent attention: Zhong et al. (2003) simulated the effects of lithospheric thickness variations on GIA in global 3D numerical models. Paulson et al. (2005) were among the first in building 3D viscosity from seismic tomography for GIA calculations. Latychev et al. (2005) investigated the impact of variations in lithospheric thickness on GIA and plate motion, followed by Klemann et al. (2008), who discussed the impact of lithospheric thickness variations and plate boundaries on the mobility of continental plates. Kendall et al. (2006) compared predicted present-day RSL rates of 3D models and a 1D model and found differences of several mm/yr. Whitehouse et al. (2006) and Steffen et al. (2006) investigated 3D and 1D GIA models in Fennoscandia and compared the results to present-day uplift rates and RSL data. Martinec and Wolf (2005) discussed the impact of 3D structure on relaxation times in Fennoscandia. Li et al. (2018) searched for a 3D viscosity model that best fits various 10.1029/2021GC009853 3 of 21 observations (paleo sea-level data, uplift rate and gravity rate) in Fennoscandia and North America. Kuchar et al. (2019) showed the importance of 3D structure for the Atlantic and Gulf coasts of North America. For eight 3D structures in Greenland, Milne et al. (2018) modeled significant differences in Holocene RSL predictions of several tens of meters in comparison to 1D models, but they also emphasize a significant spread within the 3D model ensemble. Austermann et al. (2013) discussed the impact of a subducting slab on the interpretation of sea-level data at Barbados and Klemann et al. (2007) investigated the resulting asymmetry in displacement for the South-Patagonian ice field. The prominent structural dichotomy between West and East Antarctica was investigated by Hay et al. (2017), Ivins et al. (2013), and Kaufmann et al. (2005); and van der Wal et al. (2013Wal et al. ( , 2015 investigated the influence of composite non-linear mantle viscosity due to grain size and water content in the upper mantle. Recently, Austermann et al. (2021) focused on the last interglacial and demonstrated the effect of 3D Earth structure on RSL.
One drawback of 3D models are the computational costs: less than one hour to solve the 1D problem in the spectral domain on a personal computer versus days to solve a finite-element model in the spatial domain on a computer cluster. Accordingly, the validity of 1D approximations for specific regions has been investigated in several studies: Paulson et al. (2005) made a comparison between the post-glacial rebound (PGR) derived from a spectral 1D GIA model and a finite-element 3D GIA model. Results from that study suggest that near-field observations tend to depend on the viscosity structure beneath the ice load, while the observations away from the ice sheet tend to depend on the viscosity below both the ice regions and the farfield site. A et al. (2013) applied a 3D viscosity structure, 1D averages and ICE-5G history to investigate GIA models for Antarctica and Canada. They argued that the 3D model can be better approximated using a 1D Canadian average than the 1D global average and for Antarctica, they showed that a regional 1D viscosity does not work well in reproducing the RSL observables from 3D models. Crawford et al. (2018) developed a method to quantify the sensitivity of postglacial RSL to ice history and laterally varying viscosity by applying adjoint-derived sensitivity kernels. Hartmann et al. (2020) proposed a method where they combined the response of locally adjusted Earth models to infer the sea-level response, and Bartholet et al. (2020) showed that patchy regions with very low viscosities have a minor impact on the global sea level pattern. Nevertheless, further studies are required by the GIA community, especially to extend the reasonable range of 3D GIA models as discussed at a GIA workshop in Ottawa in 2019, <https:// www.scar.org/scar-news/ serce-news/ottawa-gia-workshop>, and to estimate realistically uncertainties when determining viscosity variations from seismic tomography.
In this study, we generate 3D viscosity structures based on seismic tomography and geodynamic constraints (Steinberger, 2016). To consider uncertainties in 3D viscosity distributions, we set up an ensemble of 3D Earth structures that differ in lateral and radial viscosity distribution (Section 2). The Earth structure ensemble and three different glaciation histories ICE-5G, ICE-6G, and NAICE are applied to a 3D GIA model, called VILMA. The numerical efficiency of this code based on the spectral-finite element approach (Martinec, 2000) allows us to discuss such a model ensemble (Section 3). We investigate how variations in 3D viscosity structure influence GIA focusing on RSL predictions of globally distributed sites. Furthermore, we focus on ice center (North America, Antarctica) and peripheral regions (Central Oregon Coast, San Jorge Gulf) to investigate the effect of regionally adapted 1D structures (Section 4).

Viscosity Structure Parameterization
To investigate the effect of lateral variations in viscosity structure on GIA, we follow the reference case of Steinberger (2016) and adopt 3D Earth structures derived from tomography models SL2013SV (Schaeffer & Lebedev, 2013) above 200 km depth and the 2010 update of Grand (2002) covering greater depth, with a smooth transition. The transition is placed at 200 km depth because SL2013SV is better resolved in the uppermost mantle, showing features such as cratons more clearly. Yet its resolution decreases with depth, such that the misfit between model predictions and observations, such as for dynamic topography, becomes more extensive if the transition between models occurs at greater depth (Steinberger, 2016). Below the lithosphere, seismic velocity anomalies are converted to temperature anomalies. The conversion follows model M2 in Steinberger and Calderwood (2006). Therein the relevant references are given: In the upper mantle, −(dv s /dT)/v s ≃ 10 to 15 × 10 −5 /K meaning that a negative relative seismic velocity anomaly dv s of 1%-1.5% corresponds to a temperature T that is 100 K higher. Within the lithosphere, with thickness initially also derived from the same tomography model (Steinberger, 2016), an error function temperature profile is assumed.
Temperatures are converted to viscosities, using an Arrhenius law, that is, viscosity where H is activation enthalpy, R is the universal gas constant 8.3144 J/K/mol and r is an (adjustable) factor. This activation enthalpy factor has been introduced because a non-Newtonian flow can be closely imitated by a Newtonian flow with reduced activation enthalpy (Christensen, 1983). Activation enthalpy is also adopted from Steinberger and Calderwood (2006); in the upper mantle, ∼500-700 kJ/mol (Calderwood, 1999) based on an activation energy 525 kJ/mol from Kohlstedt and Goetze (1974) for dislocation creep in dry olivine and activation volume 10-12 cm 3 /mol. These values can be effectively reduced by multiplying with r. Below the lithosphere, we apply a radial viscosity profile (RVP) (i.e., laterally averaged viscosity at a given depth) derived from pressure and radial temperature dependence and constrained by modeling geoid, radial heat flux and the "Haskell viscosity average" of 10 21 Pa s (Haskell, 1935;Mitrovica, 1996;Steinberger & Calderwood, 2006). Fitting the geoid based on tomography models with a radial viscosity structure is generally regarded as appropriate. For example, Ghosh et al. (2010) find that the geoid calculated from tomography is hardly affected by the presence of lateral viscosity variations. More precisely, Equation 1 specifies viscosity variations due to temperature deviations from an (adiabatic) temperature profile. Since these deviations also include the temperature drop in the lithosphere, we do not explicitly define a higher lithosphere viscosity in the radial profile but adopt the radial profile only below the lithosphere. The spatial resolution of the applied Earth structures extends up to spherical harmonic degree/order (d/o) 63, corresponding to a spatial resolution of around 320 km. The applied Earth structures are described in more detail in Steinberger (2016) and Steinberger and Calderwood (2006), including the conversion from tomography to temperatures and viscosities, as well as the constraints from mineral physics and surface observations.
To reduce computing time, we limited derived viscosities lower than 10 19 Pa s to this threshold value being aware that small time-scale dynamics at low viscosity zones (e.g., West Antarctica) are thus not completely resolved. Furthermore, we consider lateral variability only above 870 km depth, due to the fact that lithosphere and upper mantle are the regions with most significant variability, as documented in, for example, Mooney et al. (1998) and Artemieva et al. (2006), respectively, and because, in general, lateral variability in mantle tomography models decreases with depth (e.g., Ritsema et al., 2004). Furthermore, the sensitivity of the GIA process to viscosity structure decreases with depth (e.g., in Fennoscandia, Steffen et al., 2006). Below 870 km depth, the considered viscosity structure is defined as the mean lateral viscosity at the specific depth. Hence, we discretize the model into 114 3D layers between surface and 870 km depth and 50 1D layers from 870 km to the Earth's core. The core itself is considered as a boundary condition resembling the response of an inviscid fluid.

Variations in 3D Viscosity Structure
For the setup of the ensemble of Earth structures, we combine six different activation enthalpy factors, r, and three different RVPs. The factor r multiplied with the activation enthalpy (see Equation 1) is varied between 0.2857 and 1.0 to mimic dislocation creep (non-Newtonian viscosity) with different strain-rate dependence of viscosity (Christensen, 1983) where r = 1 implies diffusion creep. We choose six factors, 0.2857, 0.3, 0.4, 0.5, 0.75, and 1.0; the higher the factor r, the stronger is the viscosity dependence and, so, the lateral variability (see also Steinberger & Calderwood, 2006).
We consider three different RVPs, s16, sc06, and sc06b, which are applied in order to allow for different viscosity contrasts between the upper mantle and transition zone at 410 km depth and between the transition zone and lower mantle at 660 km depth (see Figures 2 and 3). The structure s16 shows a substantial viscosity increase with depth at 410 km depth, which is generally not considered in GIA models like VM2 and VM5a . The sc06 RVP offers lower viscosities in the transition zone between 410 and 660 km depth than in the layers located above and below this zone. For the RVP of sc06b, there is no viscosity step at 410 km, which coincides with the general assumption for the depth range of a homogeneous upper mantle usually applied in GIA models.
We can distinguish 18 variations of 3D Earth structures (Figures 2 and 3), which we split according to the RVPs into three classes, Class-I, Class-II, and Class-III representing the s16, sc06, and sc06b, respectively. We label each viscosity structure by v_[r]_[RVP], for example, for v_0.4_s16 a factor r = 0.4 and RVP s16 is applied. In combination with r = 0.4 -which was chosen as reference value because Christensen (1983) showed for 2D numerical experiments that the properties of non-Newtonian flow with n = 3, which is approximately appropriate for dislocation creep, can be closely imitated by Newtonian flow with activation enthalpy reduced by a factor r = 0.3 -0.5 -the structure v_0.4_s16 corresponds to radial viscosity structure "SL + Gra3" in Steinberger (2016), v_0.4_sc06 corresponds to Model M2 and v_0.4_sc06b corresponds to Model M2b, the latter two in Steinberger and Calderwood (2006).
For purpose of visualization and statistics (Figures 1 and 2), we show the lithospheric thickness T lith and the logarithmic depth-averaged viscosity value, <log η(z)> z , for three layers, asthenosphere (a), upper mantle (um), and transition zone (tz). The "lithospheric thickness" is defined as the minimum depth below which viscosities are less than 10 23.5 Pa s. For example, for reference structure v_0.4_s16, this corresponds to a temperature of 1,087 K (100 km) and 1,153 K (200 km). Steinberger (2016) uses a higher temperature of about 1,405 K, which results in larger lithospheric thicknesses. This viscosity-based definition of the lithosphere differs from that commonly used in GIA, where the lithospheric thickness is considered as an effective elastic thickness for which, on GIA related loading times, the lithosphere behaves elastically. The asthenosphere considered here comprises the region between T lith and 225 km depth, approximately the lower boundary of a zone with reduced seismic velocities and increased attenuation, for example in PREM (Dziewonski & Anderson, 1981). However, calling this depth range "asthenosphere" may be misleading: whereas in areas of thin lithosphere it includes the depth range with the viscosity minimum, this is not the case in areas of a thick lithosphere, where the gradual viscosity decrease below the lithospheric mantle exceeds the depth of 225 km. The "upper mantle" comprises the region between 225 km (or T lith for structures with lithosphere >225 km) and 410 km depth defining the discontinuity to the mantle transition zone. The "transition zone" is specified between 410 and 670 km depth, the top boundary of the lower mantle. Note that we are using the logarithmic viscosity, hence, also the ensemble mean and ensemble standard deviation are derived from the logarithmic viscosity. Table S2 in Supporting Information S1 shows the maximum and minimum values of the ensemble mean and ensemble standard deviation. In the data publication , we provide all 18 3D Earth structures.
From Figure 1, we clearly distinguish tectonic regimes of divergent plate boundaries or young lithosphere characterized by thin lithosphere (<30 km) and low viscosities in the asthenosphere and upper mantle (∼10 20 Pa s) (e.g., Cascadia subduction zone, West Antarctica), while regions of deep rooting cratons are characterized by thick lithosphere (>200 km) and high viscosities (>10 23 Pa s) in the asthenosphere (e.g., Fennoscandia, North America, East Antarctica). In general, larger values in lithospheric thickness or viscosity cause more considerable variations between the models (larger ensemble standard deviation) (Figure 1) resulting in an ensemble variability that depends on the region (Figures 3d-3h). The rather large standard deviations in the mantle viscosities of the whole model ensemble (Figures 1f-1h) are dominated by the three chosen RVPs. This is also visible in Figure 2, where we present the lateral variability by the median and the 5th, 25th, 75th, and 95th percentile of the respective layers. The figure reveals general information about the variability depending on the activation enthalpy factor and how the RVP influences the differences between the three classes. It becomes evident that a larger activation enthalpy factor generates increased lithospheric thickness and larger variability in each layer. Furthermore, the asthenosphere variability is more extensive than the upper mantle and transition zone for all classes. When comparing the three classes, the viscosity contrast between upper mantle and transition zone results in an opposite behavior, which is due to the constraint associated with the Haskell viscosity average.

Derived 1D Viscosity Structures
In order to discuss the impact of lateral variability, we construct three sets of 1D structures based on the 3D ensemble. For the first set, we took the global mean of each ensemble member to derive the respective radial profile. From the resulting profiles in Figures 3a-3c, we find that the activation enthalpy factor mainly generates variations in lithospheric thickness and asthenosphere, whereas the specific RVP causes the expected variations between asthenosphere, upper mantle, and transition zone. In contrast, the variability inside each RVP class is negligible in the upper mantle and transition zone.
For the second set, we construct regionally adapted viscosity structures for two core regions we identified for this study, the Central Oregon Coast (#5) and San Jorge Gulf (#4) (see Figures 3d-3f). Due to their tectonic setting, we expect significant lateral heterogeneity as it was previously discussed in the context of 3D Earth structure (Clark et al., 2019;Klemann et al., 2007). The first region of interest is #5 with low viscosities in the asthenosphere due to its coincidence with the Cascadian subduction zone. The cross-section in Figure 4a shows the significant change toward the cratonic lithosphere below the former Laurentide ice sheet (LIS) (Figure 4b). The second region of interest is #4 at the Atlantic coast of southern Patagonia. Southern Patagonia is characterized by subduction of an immature part of the Antarctic plate below the South American plate. In addition, the Patagonian ice sheet (PIS) extended along the Chilean coast of Patagonia. The cross-section in Figure 4c visualizes the increase of lithosphere thickness and viscosity toward the South American plate.
The regionally adapted 1D viscosity structure should represent the local viscosity structures of the respective region. For this, we define a radius of interest inside which we determined a laterally averaged viscosity structure (marked as gray ellipses in Figures 4b and 4d). For #5, we chose a radius of 5° (∼560 km) and for #4 we chose two radii, 5° and 2.5° (∼280 km). We choose two radii for #4 to investigate the effect of the chosen radius on the mean viscosity structure and RSL predictions since the region is characterized by a sizable lateral viscosity variation on a small regional scale. In contrast to the negligible impact of the activation enthalpy factor r on the deeper viscosities in the case of the global means (Figures 3a-3c), the variations remain significant in the upper mantle and transition zone for the regionally adapted structures. As expected, the adapted structures show significantly thinner lithospheres and lower viscosities in the asthenosphere. Furthermore, the general viscosity increases with r due to Equation 1 does not hold for #5 structure and #4 structure with a 2.5° radius where we find smaller viscosities in the asthenosphere.
For the third set, we construct regionally adapted viscosity structures for North America and Antarctica. As a mask we choose the ICE-6G extension at 21 ka BP for the North American and Antarctic ice-sheet, respectively (Figures 3g and 3h). We calculated the lateral viscosity mean in these areas and used the 1D viscosity profiles as global Earth structures, respectively.

Model Setup
For the prediction of RSL during the last glacial cycle, we apply the VIscoelastic Lithosphere and MAntle model VILMA (Klemann et al., 2008;Martinec et al., 2018). In this code, the field equations for an incompressible self-gravitating viscoelastic sphere are solved in the spectral-finite element formulation of Martinec (2000). Therein, the sea-level equation after Kendall et al. (2005), including deformation and gravitational effects, rotational feedback, floating ice and moving coastlines is solved following Hagedoorn et al. (2007) and Martinec andHagedoorn (2005, 2014). Elastic parameters and the density structure of the san-jorge_rad5, (f) san-jorge_rad2.5, (g) north-america, (h) antarctica. For comparison, 1D viscosity structures VM2  and VM5a  are shown.
Earth are derived from the Preliminary Reference Earth Model (PREM, Dziewonski & Anderson, 1981). Note, that for an incompressible stable sphere, a constant density would be sufficient, but as we consider the buoyancy forces internally and at each density contrast, we have to use a realistic density structure. The field equations of the linearly viscoelastic continuum, as well as the sea-level equation, are solved by an explicit time-differencing scheme, which demands a timestep of two years due to the lower viscosity threshold of 10 19 Pa s. The spatial representation of the sea-level equation at the surface, as well as the lateral variations in the 3D layers inside the continuum, are solved on a Gauss-Legendre grid of 256 × 512 grid points (n128), which is consistent with the spectral resolution in spherical harmonics up to Legendre d/o 170, corresponding to a wavelength of ∼120 km. This resolution is by a factor of 3 larger than the resolution of the 3D viscosity structure under consideration here, which only is up to d/o of 63. However, it is necessary for the calculation of the deformation response that includes the discretization of the ice load and the solution of the sea-level equation for changing bathymetry and shorelines. The radial finite-element node-distance amounts to 5 km for the upper 420 km depth, 10 km between 420 and 670 km and 40-60 km between 670 and 6,371 km. In order to reach a present-day sea-level corresponding to the present topography, four integration loops of the whole glaciation history were performed for each ensemble member. To consider a range of different surface loading histories, we implemented three published glaciation histories (Text S1 in Supporting Information S1) covering the last glacial cycle: ICE-5G , ICE-6G ) and the regional NAICE . ICE-5 G and ICE-6G are not independent models and this is therefore not a rigorous exercise, but we nevertheless expect that this model comparison can provide some insight into the dependence of results on glaciation history.
We composed in total about 300 GIA models from different combinations of 3D and 1D Earth structures with the three ice histories (Table S1 in Supporting Information S1). ICE-5 G and ICE-6G were combined with the complete 3D ensemble, whereas for NAICE we considered only eight 3D structures ([r]_s16_3D, 0.4_sc06_3D and 0.4_sc06b_3D). We label the models with m_[r]_[RVP]_[structure]_[ice], according to the considered Earth structure and the respective ice history. The 3D model ensembles (44 members) we label as 3D_ICE-5G, 3D_ICE-6G and 3D_NAICE. For the numerically much cheaper 1D Earth structures, we considered all possible combinations, that is, 54 models for the global mean structures and 162 for the three regionally adapted structures oregon_rad5, san-jorge_rad5 and san-jorge_rad2.5, where rad5 and rad2.5 refer to the respective radius in a degree that is used to create the regionally adapted models. Furthermore, we combined the 36 1D regionally adapted structures for North America and Antarctica with ice history ICE-6G.

Relative Sea-Level Prediction
In order to analyze the deformational behavior depending on the viscosity structure, we use the RSL. We focus on the RSL since 14 ka BP, as an increasing amount of sea-level data is available to compare a significant part of the deglaciation phase with our model output.

Global Range of RSL Predictions
The RSL variability due to the 3D Earth structure is shown in Figure 5, where the ensemble range of the predicted RSL is plotted at 14 ka BP for the ensembles 3D_ICE-5G and 3D_ICE-6G. As expected, the largest variability in RSL occurs in regions covered by the dominating Laurentide ice sheet LIS, the Greenland ice  sheet GIS, the European ice sheet EIS, extending from the British Isles to the Kara Sea, and the Antarctic ice sheet AIS, but also in smaller ones like Iceland and the Patagonian ice sheet PIS. Depending on the ice history, the variability due to Earth structure variations has different emphasis (Figures 5a and 5b). For example, 3D_ICE-5G shows large variability at the Antarctic Peninsula and Ellsworth Land, while for 3D_ICE-6G, large variability occurs in Marie Byrd Land. The ensemble 3D_NAICE reveals that, due to less ice volume in North America, the resulting RSL variability is reduced in this region ( Figure S1 in Supporting Information S1). These results are comparable to Li et al. (2020) who focused on North America also considering the ice history ICE-6G, where our ensemble shows a 10% smaller mean RSL at 15 ka BP, but inside the ensemble range of 20%. Figure 6 shows the variability of RSL at 14 ka BP for the eight locations indicated in Figure 5, split into the three ice histories. Because the NAICE ice sheet only covers North America, we show the corresponding results only for #1 and #5. To visualize the range of RSL predictions compared to data uncertainty we apply paleo sea-level data (Text S2 in Supporting Information S1). The visual comparison of the sea-level indicators (SLIs) and predicted RSL curves for 3D and derived 1D global mean models between ∼15 ka BP and present day at #1, #2, #3, #4, #5, and #8 (Figures 7a-7f and Figure S3 in Supporting Information S1) shows that our models are able to reproduce reasonably RSL observations.

RSL for 3D Models and 1D Global Mean Models at 14 ka BP (Figure 6)
For Churchill at the LIS, the considered 3D model ensembles reveal a broad range of RSL with ∼150 m for 3D_ICE-5G, ∼100 m for 3D_ICE-6G and ∼50 m for 3D_NAICE. The smaller RSL height is related to the smaller total amount of this model's LIS magnitude and the smaller ensemble size of 3D_NAICE. For Angermanland at the EIS, the RSL of 3D_ICE-5G and 3D_ICE-6G are comparable with ∼150 m range. For the Ross Sea, the ranges are 26 m (3D_ICE-5G) and 6 m (3D_ICE-6G). For Oregon related to the LIS, the amount of the ranges is 23 m (3D_ICE-5G), 13 m (3D_ICE-6G), and 15 m (3D_NAICE) and, for the San Jorge Gulf related to the PIS, ranges are similar. For the far-field locations, the RSL range is smaller, but still significant with around 5 m in Rao-Gandon Area, Singapore and Pioneer Bay.
If we compare the RSL of the 3D models with that of the corresponding 1D global mean models we observe significant deviations in the ensemble mean. In detail, we observe a much larger difference for locations #2 and #3 than for #1. This might be because the ice sheets that melted in #2 and #3 were smaller, hence more sensitive to shallower depth levels, where the lateral variability of viscosity is more pronounced. Europe has a comparatively thick lithosphere and high viscosity, hence the 3D models give smaller uplift than 1D models, whereas it is the other way round in Ross Sea, where the lithosphere is thin and sublithospheric viscosity is low. On the other hand, at #5 there is less subsidence for 3D models than for 1D models, despite a thin lithosphere and a low asthenosphere viscosity. A reason might be that for the 3D model, the marginal regions are less well coupled to, and hence less affected by the ice sheets.

Influence of Activation Enthalpy Factor r and Radial Viscosity Profile RVP
The models reveal that the RSL variability is regionally and temporally variable (Figures 6 and 7). For example, in site #2 (Figure 6), viscosity variations due to different activation enthalpy factors r cause more RSL variability than changing the RVP. In site #1, variations in r cause less variability, although, for site #1, higher absolute values for RSL are predicted. The amount of variability due to the factor r may arise due to corresponding variations in the transition-zone viscosity (compare Figures S2a and S2b in Supporting Information S1).
Depending on location, the highest RSL are predicted by Class-II models and the smallest by Class-I models (#1, #2) or the other way around (#4, #5), while RSL predictions of Class-III models lie in the middle of predicted RSL at 14 ka BP ( Figure 6): For sites #1 and #2, which are characterized by a thick lithosphere, the larger RSL for Class-II models might be related to the low viscosity in the transition zone resulting in a more efficient GIA. At subsiding peripheral sites #4 and #5, it is more complex. Mostly, Class-I models predict a highest RSL (i.e., least subsidence) despite of lower asthenosphere and upper-mantle viscosities and a thinner lithosphere. This might be because in this case the locations are less well coupled to, and hence less affected by the ice sheets. We can also see that the behavior can vary with time (e.g., site #5 at 16 ka BP, Figure 7e), which might hint at a superposition of the effects from the ice center and peripheral regions.
At sites #1, #2, and #3, the increasing activation enthalpy factor r results in lower RSL predictions. The larger values of r lead to the thicker lithosphere and higher viscosity values in the asthenosphere ( Figures  S2a-S2c in Supporting Information S1), and this causes less efficient and hence smaller uplift. At sites #4 and #5, it is more complex. Depending on time, the effect of variations in factor r can result in opposite RSL predictions at one location (e.g., Figure 7e, 14-16 ka BP). The effect of r on the Earth structure changes with depth and between both locations. For example, in the transition zone and upper mantle under site #5, a larger factor leads to lower viscosities, while under site #4 the opposite is true (Figure 3 and Figure S2 in  Supporting Information S1). In the asthenosphere below site #4, higher factors lead to lower viscosities. In contrast, site #5 stands out for a substantial viscosity jump within the asthenosphere with a changing influence of r on the viscosity. We assume that the complex relationship between factor r and viscosity structure beneath the peripheral locations is not the only reason for the RSL behavior. The temporal change in the effect of factor r and RVP at the peripheral locations might indicate that also the laterally variable structure in the surrounding influences the RSL predictions.

Regionally Adapted Models for Patagonia and Oregon
We chose San Jorge Gulf (#4) and the Central Oregon Coast (#5) to test the impact of regionally adapted structures due to their tectonic setting and strong lateral heterogeneity in mantle structure. Site #4 at the southeastern coast of Argentina is located east of the PIS with ice thicknesses of several hundred meters during the last glacial cycle (Hulton et al., 2002; and recent ice mass flux during the late Holocene (Rignot et al., 2003). This region, south of the Chilean Triple Junction, can be characterized tectonically as an immature ridge of the Antarctic plate subducting below the South American plate (Bird, 2003;Mc-Culloch et al., 2000). Furthermore, the continent is characterized by a thin lithosphere and low viscosities below the overriding South American plate (Richter et al., 2016), whereas a large increase in lithospheric thickness and viscosity occurs around 500 km offshore the eastern coast (Figures 4c and 4d). An essential aspect of the GIA process is that large areas of the Atlantic shelf have been exposed during glaciated times, controlling the hydro-isostatic adjustment process in this region. The predicted RSL of the regionally adapted ensembles (san-jorge_2.5 and san-jorge_5) are plotted in Figures 7g and 7h. The impact of the radius seems to be small for the majority of Earth structures. But in comparison to san-jorge_5, Class-I model with the highest activation enthalpy factor predicts a higher RSL for san-jorge_2.5 (∼5 m at 8 ka BP), which might be associated with lower viscosities and less smooth viscosity changes in the asthenosphere (Figure 3). However, both 1D regionally adapted ensembles as well as the 1D global mean ensemble for both ice histories, underestimate the 3D RSL predictions which shows strong evidence for 3D structure in this region. The surrounding rebound may be affected by the viscosity dependent response of the solid Earth due to the PIS and the AIS.
In correspondence to site #4, the site #5 is formed by the Cascadian subduction zone where the Juan de Fuca, Gorda, and Explorer plates subduct below the North American Plate (Schmandt & Humphreys, 2010). The high-viscosity material of the thick North American continental lithosphere in the eastern part borders on the low-viscosity material of the oceanic lithosphere in the western part, resulting in a large viscosity gradient (see Figures 4a and 4b). The Cordilleran ice sheet did not cover site #5 but influenced this region in concert with the LIS. The model results of RSL show that the thin crust and the weak mantle in this region have a strong impact on the sea-level change. Although the regionally adapted 1D models (Figure 7i and Figure S4c in Supporting Information S1) as well as the 3D_ICE-5G and 3D_ICE-6G models stay reasonably close to data points, the 1D models cannot reproduce the 3D model predictions.
have to be increased to fit sea-level data when following the central Pacific coast several hundred km to the south. Also Paulson et al. (2005) showed that 1D GIA models cannot reproduce 3D models in distal regions, because the viscoelastic response in regions that are not located beneath the ice load depends on the viscosity structure beneath the ice load as well as on the local structure. Paulson et al. (2005) showed that the RSL predictions in North America are mostly sensitive to regional 1D viscosity below North America. Motivated by these results, we extended our study by constructing regionally adapted models for North America and Antarctica to investigate if any inferred 1D viscosity structure reflects some horizontally averaged properties of the 3D mantle below the ice. In Figures 7j-7l we show the RSL prediction for the locations Churchill, Oregon, and Ross Sea for regionally adapted 1D models for the North American ice-sheet and the Antarctic ice-sheet. In agreement with Paulson et al. (2005) the RSL predictions from the regionally adapted models tend to reproduce the 3D model predictions at Churchill (Figure 7j). The relaxation behavior is similar between the models, but in detail, depending on the Earth structure, 3D and 1D models can reach significant differences in magnitude. For example, the predictions from models m_0.2857_sc06 are very similar. The predictions from models m_1.0_s16 differ by about 15 m at 11 ka PB. The differences between 3D and 1D model predictions might be due to the fact that also within the averaged North American Earth structure lateral variations occur and that the regionally averaged 1D structure ( Figure 3g) differs from the radial structure directly beneath Churchill of the 3D structure ( Figure  S2a in Supporting Information S1). A higher activation enthalpy factor correlates with higher differences between 3D and 1D model RSL predictions (Figure 7j). Structures with higher factors are associated with larger lateral variations (Figure 2), which would support the assumption that the lateral variations within the North American craton influence the RSL predictions.

Regionally Adapted Models for North America and Antarctica
For Oregon (Figure 7k) located in the peripheral forebulge of the North American ice-sheet (Cordilleran part), the 1D and 3D models differ significantly, which is due to the strong lateral contrast in the Earth structure, and which was already shown by the Oregon regionally adapted model (Figure 7i). This supports the results by Paulson et al. (2005). Both regionally adapted 1D models cannot reproduce the RSL from the 3D models indicating that the RSL is sensitive to both North American structure and the regional structure below the (peripheral/far-field) site.
Furthermore, in this region, the influence of the activation enthalpy factor on RSL predictions is vice versa for the 3D and 1D models ( Figure 7k): While for the Earth structure beneath the ice-sheet (Figure 3g and Figure S2a in Supporting Information S1), a large factor correlates with high viscosities in the asthenosphere and upper mantle, for Oregon a large factor correlates with low viscosities in the upper asthenosphere and upper mantle. The RSL predictions are influenced by a fast rebound in Oregon and a slow rebound from the North American craton. The RSL predictions from the 1D ice region models are associated with a too slow rebound by high viscosities and a thick and stiff lithosphere. For Antarctica (Figure 7l), all 1D models underestimate the 3D predictions. Due to the very strong contrast in Earth structure between West and East Antarctica no 1D model can reproduce the 3D RSL predictions, which is also supported by A et al. (2013).

Discussion
In summary, our results show that depending on ice history, time, location and viscosity (activation enthalpy factor, radial viscosity profile), the RSL predictions can be influenced, which is consistent with previous studies. The non-intuitive, complex behavior between viscosity changes and RSL predictions has also been shown in Crawford et al. (2018). They calculated sensitivity kernels on defined locations (e.g., Canada, Tahiti), to evaluate the effect of viscosity changes as a function of depth and time on RSL predictions.
One interesting aspect is that the ice history NAICE improves the fit of 1D models ( Figure S3h in Supporting Information S1). In the construction of NAICE, observation data from Oregon  were considered, and the ice history was modified accordingly to fit this data. In principle, this might be a point against the demand for 3D models in GIA. For such an ill-posed problem, with many degrees of freedom concerning ice-sheet distribution, it is likely that a reasonable glaciation history can be validated with an appropriate set of sea-level data. On the other hand, we could reach such a fit without changing the ice history but considering a viscosity reduction in this region based on independent geodynamic constraints. From the observational point of view, one cannot favor one approach over the other.
Furthermore, viscosity depends on the forcing timescale. Therefore, constraints from geodynamics might lead to slightly different viscosity structures than is appropriate for GIA (Ivins et al., 2020;Lau & Holtzman, 2019). The consideration of a more complex rheology may be important when bridging effective viscosities from GIA to post-seismic relaxation or seismic wave attenuation. However, the scaling of our viscosity structure is based on the Haskell viscosity average of 10 21 Pa s, already derived from GIA. So, the application of such a rheological model might be a good starting point when considering a seismic tomography model which also contains an attenuation model to constrain the transfer function from seismic-velocity to viscosity variations (Benjamin et al., 2006).
The viscosity structures of this study are constrained by the geoid. Note that even a 1D viscosity model for the earth remains poorly constrained by the observations of the geoid or GIA, especially for multiple layers of viscosity. Paulson et al. (2007) showed that it is difficult to use GIA observations to constrain more than two layers of viscosity (i.e., upper vs. lower mantles). Similar studies were done for the geoid; Thoraval and Richards (1997) also found it hard to constrain details of viscosity structure. Besides the three radial viscosity profiles used in our study, other possible viscosity structures would fit the geoid (Forte & Mitrovica, 2001;Hager & Richards, 1989;Liu & Zhong, 2016;Rudolph et al., 2015Rudolph et al., , 2020. We presented two viscosity variation types to estimate the RSL variability. Other choices can produce viscosity uncertainties in 3D GIA modeling. For example, the choice of seismic tomography model controls the 3D velocity structure. While most recent tomography models agree on the main tectonic features like high velocity cratons and low velocities at mid-ocean ridges, they differ in many details, where models based on surface waves generally capture more detail in the upper ∼400 km. This leads to differences in the inferred viscosity structure and lithosphere thickness pattern as was discussed, for example, by Steinberger and Becker (2018), and can strongly influence GIA models.
Our results and data may support future work exploring the 3D Earth structures and can be used as input for further GIA models, for example, to identify regions of poorer fit to improve ice models. Furthermore, the ensemble can be used for global or different regional misfit estimations under consideration of further observations, for example, aiming to constrain the Earth structure ensemble.

Conclusions
To investigate the impact of lateral variability in mantle viscosity on RSL predictions, we considered 3D viscosity structures derived from seismic-tomography models and constrained by the geoid, the heat flux and mineral physics data from which a temperature distribution was derived. For the conversion from temperature to viscosity, a global mean viscosity of 10 21 Pa s ("Haskell viscosity average") was assumed. Different mean radial viscosity profiles (RVP) were considered, which resulted in viscosity varying by a factor of four and lithospheric thickness varying by tens of km. Lateral viscosity variability due to temperature changes was controlled by the activation enthalpy factor resulting in variations between two and four orders in magnitude. The resulting 3D viscosity structures served as parameterization of the numerical GIA code VILMA, which allows modeling the effect of lateral heterogeneity in viscosity on the flow and deformation patterns at the surface and in the Earth interior. The code considers GIA on a self-gravitating, incompressible and viscoelastic continuum in the spherical domain as well as solving the sea-level equation given rotational perturbations, moving coastlines and floating ice. As a forcing, three published ice histories, the global ICE-5G and ICE-6G, as well as the North American NAICE (embedded in ICE-6G) were considered. The resulting RSL predictions were discussed in general and, at eight specific sites that sample the spatial and temporal variability of sea-level change due to GIA.
The RSL prediction variability at specific sites, Churchill and Angermanland (central ice-sheet region), Ross Sea (near-field), San Jorge Gulf and Central Oregon Coast (peripheral region), Rao-Gandon Area, Singapore and Pioneer Bay (far-field) shows a complex dependence on activation enthalpy factor and RVP (viscosity structure) which is also influenced by time, location and glaciation history. At ice center locations we can observe that Class-I models (characterized by lower viscosities in the upper mantle than in the transition zone) and a higher factor r generally predict lower RSL, while at peripheral locations the behavior varies with time. The widespread of RSL projections in comparison to the uncertainty in sea-level data as well as the systematics between parameterization and RSL response may indicate the need for further constraints in 3D GIA modeling.
Regionally adapted 1D models for the ice center region North America have shown that they tend to reproduce the RSL predictions of the 3D models, but in detail, some ensemble members can differ significantly. Regionally adapted 1D models for the peripheral locations Central Oregon Coast and the San Jorge Gulf as well as for the ice center region Antarctica have shown that they cannot reproduce the RSL predictions of the 3D models, which is also supported by Paulson et al. (2005) and A et al. (2013).
GIA is not a local phenomenon as lateral mass transport inside the Earth mantle results in uplift and subsidence. This is especially the case for large ice sheets like the Laurentide one, where the mass transport can be described as a channel flow in the less viscous layer of the upper mantle (e.g., James & Morgan, 1990;Klemann et al., 2008) where the flow pattern is governed by the viscosity distribution. The resulting deformation at the surface is a superposition of the individual solid Earth response in the depth layers. The RSL prediction and GIA response at a location is therefore always a superposition of the local and regional responses. A regionally adapted 1D model representing the Earth's structure at the Central Oregon Coast may reproduce the localized hydro-isostatic adjustment at the Oregon Coast, as well as the near Cordilleran icesheet changes, but it will fail to reproduce the response due to the Laurentide ice sheet. Furthermore, lateral variability is justified by seismological as well as geodynamic constraints, and we note that the adapted 3D models, which were not constrained by GIA information other than the Haskell viscosity average, provide reasonable fits to the data.

Data Availability Statement
Numerical simulations were performed on the Mistral supercomputer from the German High-Performance Computing Centre for Climate and Earth System Research (DKRZ) in Hamburg. The figures are created using the GMT graphics package (Wessel & Smith, 1995, 1998. Supporting Information S1 attached to this manuscript provides Text S1 and S2, Figures S1-S4, and Tables S1 and S2. The 3D Earth structures are available at , and the predicted RSL and sea-level data are accessible at .