The Impact of a 3‐D Earth Structure on Glacial Isostatic Adjustment in Southeast Alaska Following the Little Ice Age

In Southeast Alaska, extreme uplift rates are primarily caused by glacial isostatic adjustment (GIA), as a result of ice thickness changes from the Little Ice Age to the present combined with a low‐viscosity asthenosphere. Previous GIA models adopted a 1‐D Earth structure. However, the actual Earth structure is likely more complex due to the long history of subduction and tectonism and the transition from a continental to an oceanic plate. Seismic evidence indeed shows a laterally heterogenous Earth structure. In this study, a numeral model is constructed for Southeast Alaska, which allows for the inclusion of lateral viscosity variations. The viscosity follows from scaling relationships between seismic velocity anomalies and viscosity variations. We use this scaling relationship to constrain the thermal effect on seismic variations and investigate the importance of lateral viscosity variations. We find that a thermal contribution to seismic anomalies of 10% is required to explain the GIA observations. This implies that non‐thermal effects control seismic anomaly variations in the shallow upper mantle. Due to the regional geologic history, it is likely that hydration of the mantle impact both viscosity and seismic velocity. The best‐fit model has a background viscosity of 5.0 × 1019 Pa‐s, and viscosities at ∼80 km depth range from 1.8 × 1019 to 4.5 × 1019 Pa‐s. A 1‐D averaged version of the 3‐D model performed slightly better, however, the two models were statistically equivalent within a 2σ measurement uncertainty. Thus, lateral viscosity variations do not contribute significantly to the uplift rates measured with the current accuracy and distribution of sites.

. Average GPS uplift rates (mm/yr) over the periods (a) 1992-2003 and (b) 2003-2012 derived by . Two uplift peaks with rates >30 mm/yr can be seen at the Yakutat and Glacier Bay ice fields. Mapped quaternary faults are obtained from the Alaska Quaternary fault and folds database (Koehler, 2013).   Elliott et al. (2010) and subsequent studies increased the spherical harmonics up to degree and order 2,048. This way, small ice load changes and their effects could be resolved with higher spatial detail. In addition, Elliott et al. (2010) used a refined GPS data set with higher accuracy and density. This mainly resulted in constraints with a thinner lithosphere. Sato et al. (2011) investigated increased PDIM rates using the icerate model by Larsen et al. (2007), which resulted in a larger asthenospheric thickness being required to fit the data.  modeled higher PDIM rates, derived from the ice-rate model from Berthier et al. (2010). Overall, a thin lithosphere (50-70 km) underlain with a low-viscosity asthenosphere (2.5 × 10 18 -3 × 10 19 Pa-s) is preferred. However, some areas remain either underpredicted (e.g., the YK ice fields) or overpredicted (e.g., Haines to Juneau) . These discrepancies are likely due to systematic errors in the ice load model or the Earth model.
The 1-D parameterization in previous GIA studies may not accurately represent the true 3-D Earth structure in this region which may affect the GIA predictions. Studies of other regions have shown that a 3-D structure has a large impact on the GIA predictions (e.g., Li et al., 2020;Spada et al., 2006;. Global seismic tomography studies (e.g., Schaeffer & Lebedev, 2013) and a regional tomography study for Alaska by Jiang et al. (2018) show that lateral variations in seismic velocities exist in the region, which likely correspond to variations in temperature and composition. The actual 3-D structure in Southeast Alaska is more complex due to the geologic history of tectonism. Before 40-60 million years ago, depending on the plate reconstruction considered, this region was a subduction zone (Engebretson et al., 1984;Fuston & Wu, 2020;Haeussler et al., 2003). Since that time, the region has been subjected to shear deformation and substantial margin-parallel (northward) transport of material via strike-slip faulting (DeMets & Merkouriev, 2016). Offshore of Southeast Alaska, there is an abrupt transition from continental to oceanic lithosphere across the Queen Charlotte fault.
Seismic anomalies can be converted to viscosity variations. For example, Ivins and Sammis (1995) used a scaling relationship for the conversion, by relating density anomalies to temperature anomalies based on the notion that seismic anomalies are caused by temperature variations alone. In reality, non-thermal effects such as compositional heterogeneity can also affect seismic anomalies. Non-thermal effects can also play a role in continental regions characterized by iron depletion or in tectonically active regions (i.e., characterized by partial melt and/or high water content) (Artemieva et al., 2004). In the GIA studies by Wang et al. (2008) and Wu et al. (2013), the scaling relationship of Ivins and Sammis (1995) is multiplied by a scaling factor, which represents the fractional thermal contribution to seismic anomalies. Wu et al. (2013) found a best fit to GIA observations for the case that thermal effects have a dominant control over seismic anomalies beneath Fennoscandia, with 65% control in the upper mantle which increases with depth into the lower mantle. However, uncertainties related to the relative thermal contributions increase with depth because the GIA process in Scandinavia is mostly sensitive to the upper mantle. Other methods to determine viscosity rely on flow laws and an estimate of mantle temperature (e.g., . In this work, we will mainly use the method of Wu et al. (2013) to constrain the thermal effect on seismic anomalies by using seismic tomography and GIA due to post-LIA ice load changes. Because of the short wavelength of the ice load, we restrict the investigation to parameters in the shallow upper mantle (<400 km), since crustal velocities are not very sensitive to viscoelastic relaxation in deeper mantle layers . A secondary aim of this research is to determine the importance of 3-D viscosity structure for this regional loading problem, considering the high computational cost of 3-D models. The scaling factor, which represents the fractional thermal contribution, is also a measure for the magnitude of viscosity variations. Therefore, the obtained scaling value will reveal to what extent 3-D viscosities play a role in the area. Predictions for 3-D models are compared to those of 1-D models.
To summarize, we aim to answer the following research questions: 1. How do GIA models constrain the regional 3-D viscosity structure? 2. What is the thermal contribution to shear wave anomalies in the region? 3. What is the effect of lateral viscosity variations on GIA model predictions?
The GIA modeling was performed with a Finite Element (FE) model which allows lateral heterogeneity. Section 2.1 of this paper explains the FE model setup and the Earth model parameters. Section 2.2 briefly describes the ice load model and Section 2.3 describes the method to retrieve viscosities from seismic 10.1029/2021JB022312 4 of 17 tomography. The suite of viscosity models explored in this area are presented in Section 3.1. In Section 3.2, the model misfits are evaluated using a chi-square test. The role of 3-D viscosity variations is evaluated in Section 3.3. An alternative approach to retrieving a 3-D viscosity distribution is throughflow laws for olivine, where 3-D variations result from variations in temperature. We evaluate this approach in Section 3.4. To compare our results with earlier studies, which are all based on incompressible earth models, we address the role of material compressibility in Section 3.5. Implications of our results for the regional Earth structure are discussed in Section 3.6. Model limitations are described in Section 3.7, followed by the conclusions.

Model Setup and Geometry
In this research, the FE method is used to model deformation and stress in the Earth as implemented in the commercial FE package ABAQUS FEA (Hibbitt et al., 2016), following the approach by Wu (2004). The GIA model in this research adopts a flat-Earth approximation. The validity of the flat-Earth approximation was shown in Wu and Johnston (1998) for loads up to the size of the Fennoscandian ice sheet. Hence, the flat-Earth assumption is reasonable considering the smaller extent of the ice load in Alaska since the LIA. In addition, material compressibility is assumed, while the density perturbations associated with volumetric deformation are neglected. Compressibility effects on the buoyancy force are expected to mainly affect the horizontal displacement (Tanaka et al., 2011) which is not considered in this study. Moreover, self-gravitation is neglected. Amelung and Wolf (1994) showed that the effect of neglecting self-gravitation is partly counteracted by the flat-Earth approximation, which was also confirmed in Spada et al. (2011) and Schotman et al. (2008).
The incompressible flat-Earth FE model has been benchmarked against the normal-mode (NM) model of Hu and Freymueller (2019) (see Supporting Information S1). The FE and NM models show good agreement, where most of the difference are smaller than 1 mm/yr. The largest differences (up to 2.5 mm/yr) are near the YK ice fields, which are likely due to smoothing of the ice load model in the FE grid and approximation in the FE models. In addition, tests were performed on the FE model resolution with 10 and 15 km. The NM model uses spherical harmonics with maximum order and degree of 2,048 (∼10 km). The 10 km FE resolution model did not yield significantly better results than the 15 km resolution test (differences less than 0.5 mm/yr) and resulted in much longer computational times. For that reason, the 15 km resolution was used in further simulations as it was adequate to represent the observed deformation. For further details on the model the reader is referred to the Supporting Information S1.
The model geometry is based on work by Schotman et al. (2009) and Barnhoorn et al. (2011). The loading area consists of 155 × 95 elements with the above-mentioned resolution of 15 × 15 km. Deeper layers (starting from 90 km) have a coarser resolution: 31 × 19 elements of 75 × 75 km. Figure 2 shows the model surface geometry. The total surface area of the model is 20,000 × 20,000 km and the model extends to a depth of 10,000 km in order to minimize boundary effects. In total, 26 FE layers are created, which gives a total of 198,530 elements. The bottom and vertical edges are prescribed with boundary conditions such that the bottom edge is fixed, and the sides are limited to vertical translation. Winkler foundations (Wu, 2004) are applied at the Earth's surface and internal boundaries where density jumps occur to simulate buoyancy forces.
The shallower layers of the model have a higher resolution, whereas the deeper layers have lower resolution. The resolution of these two parts is chosen such that an even number of elements of the higher resolution part fit exactly on the top surface of the lower resolution part. Considering the mismatch in nodes and the fact that ABAQUS does not provide a standard element to model this, tie constraints are applied at the two surfaces where resolution jumps occur. The tie constraints allow for all active degrees of freedom to be equal for both surfaces. The two surfaces are defined by the upper and lower element surfaces of the two layers, respectively. The outer vertical edge elements are not taken into account as these element nodes already have a fixed constraint.
The Earth model parameters are described in Table 2. The density and shear modulus are derived from volume-averaging the PREM model (Dziewonski & Anderson, 1981 (Pasyanos et al., 2014) to infer variations in the density and elastic structure and showed that lateral variations in these parameters have a small effect on the elastic uplift. We therefore do not include a laterally heterogeneous density and elastic structure. Poisson's ratio varies with depth between 0.26 and 0.30.
The upper three layers (depth <40 km) are considered fully elastic. Below 400 km depth, the viscosity is only varying with depth. Sato et al. (2011) showed that GIA is less sensitive to viscosity values below this depth due to the short wavelength of the ice load involved. Therefore, we consider it reasonable to use viscosities from a global reference model below 400 km depth and the VM5a model (Peltier et al., 2015) was adopted for this purpose. Below 40 km and above 400 km depth, each individual element within a layer is assigned to an individual viscosity value.

Ice Load Model
Upon ice removal, the Earth responds with an instantaneous elastic response and a time-delayed viscous response. The timescale associated with the viscous flow is related to the characteristic relaxation time of the mantle. Due to the presence of low viscosities in the asthenosphere, the associated relaxation times are decades to years, which are comparable to the timescales of recent ice loss and make it difficult to separate elastic and viscous processes. For that reason, both LIA and PDIM glacier variations are modeled together in the ice model. The ice load model is adopted from . We shortly summarize its main characteristics here, but the reader is referred to  for further details. The ice load model is assumed to be a function of space and time, where the ice evolution spans the last 2 ka. In essence, the ice load model consists of three sub-models defined for selected areas: the whole region, GB and the YK ice fields. The last two are necessary because ice mass loss was asynchronous with respect to the regional model. The GB ice field experienced a large ice volume loss about 200 years ago and the YK ice field experienced an accelerated ice mass loss during the last two decades. The reader is referred to Figure 5b in  for the time history of the ice loads. The glacier evolution is essentially the same as in , except for the adoption of the late twentieth century ice rate map from Berthier et al. (2010) (here referred to as the Berthier model). The data used for the Berthier model covers the period between 1962 and 2006. In their ice model,  assume these rates to represent the average ice wastages between 1900 and 1,995. For PDIM rates (1,995-2,012) the Berthier model is extrapolated and ice wastage is enhanced by a factor of 1.8 and 2.2 for the periods 1995-2003 and 2,003-2,012, respectively. PDIM has a contribution to the current uplift rates of up to 45% and 25% for YK and GB, respectively (see Figure S7 in Supporting Information S1), which is larger than the values found in . This is largely attributed to the higher PDIM rates modeled.
The regional ice load model is given by 677 disks with radii between 10 and 11 km, whereas the GB ice load model is represented by 5 additional disks with radii between 13 and 19.5 km. The disks around YK, which are given by the regional model, are subjected to ice loss rates three times larger than the regional model. This loading is interpolated to the FE grid, while conserving total mass change at each time step. The fraction of each disk covered by a rectangular FE grid element is assigned to a rectangular grid cell. The ice load distribution is shown in Figure S3 in Supporting Information S1. The mesh size needs to be small enough so that errors around the ice load edges are minimized. Our benchmark analysis (see Supporting Information S1) showed that a resolution of 15 km was sufficient to this end.

3-D Viscosity Calculations From Seismic Velocity Anomalies
The 3-D viscosity structure is estimated directly through seismic tomography as described in the approach by Wu et al. (2013). In the crust seismic anomalies are mainly controlled by composition, whereas in the upper mantle seismic wave anomalies are for a large portion controlled by temperature (e.g., Cammarano et al., 2003;. Assuming that temperature variations are responsible for viscosity anomalies, Ivins and Sammis (1995) introduced a scaling relationship between seismic velocity anomaly and viscosity that computes viscosity anomalies based on the effect of temperature and mineral physics. Wu et al. (2013) slightly modified this relationship and included a parameter to scale the viscosity anomaly based on the thermal contribution. The scaling relationship is given by : is the fractional shear wave anomaly computed with respect to the reference seismic anomaly profile  ,0 s E , and is the velocity derivative with respect to temperature accounting for both anharmonic and anelastic effects. Thus 3-D temperature variations are determined from the seismic velocity variations (scaled by the parameter  E ), and viscosity variations are inferred from the temperature variations. The absolute viscosity is then related to the background viscosity and the viscosity anomalies with    The seismic anomalies are taken from the global shear wave velocity model SL2013sv (Schaeffer & Lebedev, 2013). This model defines lateral variations in velocity with respect to a 1-D velocity profile for the mantle (<1,000 km). Although uncertainties associated with the input velocity model can influence the results (e.g., Yousefi et al., 2021), we do not consider them. Shear wave velocity anomalies in Southeast Alaska are dominantly negative (see Figure 4, upper panels), which result in lower viscosities and thus a weakening effect on the upper mantle rheology. The velocity derivatives in Equation 1 are taken from Table 20.2 in Karato (2008), which represent global averages and may introduce a bias for Southeast Alaska; anelasticity is expected to play a larger role in this area due to the higher temperatures involved . Also, if indeed the mantle is substantially hydrated, the increased water content will enhance the anelasticity effects ). If anelastic effects are not taken into account (or not enough), temperatures could be overestimated and in turn result in lower viscosities. Uncertainties related to the effect of the composition, water content and partial melt, are not considered here and may also affect the  E parameter.
The 1-D background temperature profile ( 0 E T ) used in Equation 1 is a composite of the globally averaged profile by Stacey and Davis (2008) and the regional study by . Temperatures between 0 and 39 km depth are taken from Stacey and Davis (2008) Sato et al., 2011). Therefore, the globally averaged profile is not considered suitable.
The search space consists of the background viscosity profile and the  E parameter in the upper mantle. We limit the search grid to the upper mantle as GIA does not constrain deeper viscosities due to the short wavelength of the regional deglaciation (e.g., Sato et al., 2011). Choosing the background viscosity should be done carefully. The VM5a viscosity structure is not suitable for this area due to its relatively high viscosity in the upper mantle: 10.0 × 10 21 (between 60 and 100 km) and 0.5 × 10 21 (between 100 and 420 km) Pa-s, as opposed to the range (2.5 × 10 18 -3.0 × 10 19 Pa-s) for the asthenosphere found in regional GIA studies (e.g., Sato et al., 2011). We selected the best-fit earth model of  as the baseline for comparison with our analysis: the elastic lithosphere is 55-km thick and the asthenosphere is 230-km thick with a viscosity of 3.0 × 10 19 Pa-s.

Viscosity Variations and 1-D Averaged Profiles
The mantle viscosity throughout the volume was computed following Equation 1. Lateral viscosities at selected depths for the best fit model (described in Sections 3.2 and 3.3) are depicted in Figure 4 (lower panels). The seismic anomalies are directly related to the viscosities; negative shear wave anomalies indicate larger temperatures and thus lower viscosities. The largest variations are seen at approximately 80 km and variations decrease with depth. In deeper layers, the viscosity variations decrease as seismic velocity variations decrease, because we use a constant scaling relation throughout the mantle. In addition, we ran a number of models using an alternative approach (described in Section 3.4), which entails inferring viscosities from temperature variations from the global temperature model WINTERC-G  throughflow laws for olivine . . These profiles illustrate a tradeoff between these two parameter values for models providing a good fit to the GPS observations, such that the scaling factor  E needs to be larger for models with a higher background viscosity in order to fit the data well. All of these models have relatively similar 1-D average viscosity profiles, with the lowest viscosities always located at shallow depth in the southeastern part of the model domain (Figure 4, left panels). However, a lower background viscosity combined with small lateral viscosity variations results in the best fit. The acceptable combinations for  0 E and  E are seen in Figure 6a, which are an approximation of a 2σ uncertainty based on a 20% misfit increase.  where round symbols indicate results from incompressible models. The best fit model (χ 2 = 13.7) is obtained for the scaled seismic anomalies approach with scaling factor 0.1 and background viscosity 5.0 × 10 19 Pa-s. The green outline in panel (a) corresponds to models with an increase in χ 2 < 20% with respect to the best fit model, which is an approximation of a 2σ uncertainty.
For the profiles shown in Figure 5b, the averaged viscosity is either too low or too high, resulting in poor prediction of the uplift rates. For example, the green line in Figure 5b has, on average, a lower viscosity than the best fit model, and the blue line has higher averaged viscosities. All of the viscosity models based on the WINTERC-G temperature model performed poorly (Section 3.4). The best-fit profile from this approach is indicated in Figure 5b (purple lines). This approach did not yield good fits because the thick lithosphere imposed by the global temperature model WINTERC-G is incompatible with the thin lithosphere required in this region (Figure 3). Further details on the results of this approach are discussed in Section 3.4.

Model Performance
The GIA model performance is tested against the GPS rates from Hu and Freymueller (2019), shown in Figure 1. The vertical uplift is due to the combined effect of GIA due to post-LIA, PDIM, and Pleistocene glaciations. In this study, we only modeled post-LIA and PDIM effects. We superimpose our modeled results with the effects of Pleistocene glaciations modeled in . The impacts of Pleistocene glaciations include contributions of the Laurentide ice sheet (ICE-3G), glaciers in southern Alaska (Wheeler, 2013) and glaciations of southern British Columbia and Cascadia (James et al., 2009), and are very small in this region (on the order of 1 mm/yr). These Pleistocene glaciations effects were computed by Hu and Freymueller (2019) using a NM GIA model, which may not represent the effects with a FE model with different Earth parameters. However, the effects by Pleistocene glaciations are sufficiently small that the overall difference is assumed to be negligible. We do not consider tectonic effects as these effects are expected to be small because of the largely strike-slip tectonics (Elliott et al., 2010) (e.g., <2 mm/yr for the fault system in southern California (Smith-Konter et al., 2014)). The misfits between the observed and predicted GIA rates are evaluated using the χ 2 test: where E N is the number of observations, i E o is the observed GPS rate, i E p is the predicted uplift rate (including Pleistocene glaciations, LIA, and PDIM effects) and  i E the GPS error. GPS observations are available for two periods: 1992-2003 and 2003-2012. Observations of both periods are combined to compute the best-fit value. Note that we do not include the error in the Pleistocene glaciation, LIA and PDIM loading models as there is no good error information available for these models.
The viscosities of the best fit model vary between 1.8 × 10 19 and 4.5 × 10 19 Pa-s. Viscosity maps at selected depths are shown in Figure 4 (lower panels). The 3-D viscosities of the best fit model are in the same range as found in previous 1-D GIA modeling studies (e.g., Sato et al., 2011). However, a large portion of the 3-D viscosities is larger than the previous 1-D GIA model, which is attributed to the inclusion of compressibility in our model (see Section 3.5). Figure 6a shows that incompressible models prefer a lower background viscosity than compressible models with the same scaling parameter. Present-day velocities increase for a compressible model in comparison to an incompressible model with the same earth model parameters. Thus, an incompressible model would require lower viscosities to achieve similar uplift rates as compressible model, given the loading model. Tanaka et al. (2011) showed that incompressible models underestimate the viscosity, which agrees with our results. Next, we investigate if the lateral variations for our best fit model are significant enough to be differentiated from a 1-D model with the GPS uplift rate data.

Role of 3-D Viscosities in the Upper Mantle
The optimal background viscosity and the scaling parameter are η 0 = 5.0 × 10 19 Pa-s and β = 0.1 (Figure 6a), respectively. The best-fit model has a fit value of χ 2 = 13.7 with residual mean of 0.3 ± 3.92 mm/yr. The residuals of the 3-D model are shown in Figures 7a and 7d. This 3-D model was able to reduce the residuals with respect to the best-fit 1-D model in ; most reductions are located inland and along the Lynn Canal, the inlet from Haines to Juneau. Because the 1-D model in  is incompressible, we compare the 3-D model with a "1-D averaged" model, where the viscosity structure is derived from the best fit 3-D model by averaging the viscosity in each layer over a certain area.
The area is confined to grid cells with prediction rates larger than 15 mm/yr. The viscosity profile can be seen in Figure 5a, which shows that the 1-D viscosity profile is closer to the upper bound of the viscosity variations in the 3-D structure. The 1-D averaged model has a fit value of χ 2 = 12.1. We also tested thresholds of 5, 10 and 20 mm/yr, which resulted in a change of +0.2, +0.3, and 0.0 in the χ 2 value, respectively. Thus, we conclude that the averaging is not very sensitive to the considered test areas.
The differences between the 3-D and 1-D averaged model uplift prediction rates are shown in Figures 7b  and 7e. The largest differences in the uplift occur in the southeastern corner (up to −2 mm/yr), which is where the largest variations in seismic velocity and thus viscosity are located. However, there are few observations available where the largest differences occur. We do see that the 1-D model systematically performs better than the 3-D model along Lynn Canal, the inlet from Haines to Juneau ( Figure S5 in Supporting Information S1). For the 3-D case, this corresponds to the zone where seismic velocity anomalies were most negative (where there are GPS measurements) and hence the viscosity becomes lower than its surroundings at shallow depths (Figure 4). This results in too large uplift rates for the 3-D model, while the 1-D model predictions are smaller and closer to the observations. The 3-D model systematically performs better around GB; however, the difference is within the GPS uncertainties. It is possible that, errors in the seismic velocity model make the 3-D model fit slightly worse than the 1-D model, or the contribution of temperature variation to viscosity might not be uniform over the area. Overall, most of the differences are smaller than 2σ, indicating that the 3-D and 1-D models cannot be differentiated from each other (Figures 7b and 7e). This implies that lateral variations do not impact the predicted uplift rate significantly.

Comparisons to 3-D Viscosities Derived From Olivine Flow Laws
An alternative approach to determine the 3-D viscosity variations is by using flow laws for creep of olivine. Here, we use the same methods in van der  for diffusion creep. For this approach, viscosities are inferred from lateral variations in temperature, water content and grain size are varied and a best fit with the GPS data is sought. We use the global temperature model WINTERC-G by , for which the average temperature profile beneath Southeast Alaska is shown in Figure 3. Water content and grain size shift the viscosity profile (the purple line in Figure 5b), whereas the temperature model determines the shape of the profile and the magnitude of the viscosity variations. For details on the method and considered parameter values, the reader is referred to the Supporting Information S1. Note that the WINTERC-G provides absolute temperatures; it is therefore not possible to use composite background temperatures shown in Figure 3 for the 3-D temperature field. We could have opted to use the 3-D temperatures derived from the background geotherm and seismic anomalies, but such a model would not be independent. Furthermore, the results show that the lithosphere of the 3-D models is likely too thick. Since WINTERC-G already has the largest temperatures ( Figure 3); replacing WINTERC-G based temperatures with the seismically derived temperatures would exacerbate this problem.
The fit results (χ 2 ) are seen in Figure 6b. The plot shows that similar fits can be obtained with certain combinations of water content and grain size: if we increase the water content (a weakening effect), then we need to increase the grain size (a strengthening effect). The best fit parameters of 8 mm grain size and 400 ppm H 2 O result in a fit value of χ 2 = 20.1. The best fit value of 20.1 is, however, more than 50% larger than the values obtained with the scaling of shear wave tomography approach. The uplift pattern (not shown) indicates that the predictions are diminished in amplitude and the largest residuals are seen at the uplift peaks. The larger misfits in this model are due to the thicker elastic lithosphere, which was induced by the WINTERC-G temperatures (Figure 3) when translating to viscosity. As shown in Figures 5b and 5a, sharp change in viscosity occurs at a depth of ∼100 km, whereas the other models show a thinner lithosphere in which this viscosity contrast occurs at 55 km. This approach results in larger lateral viscosity variations, where viscosity contrasts are up to a factor of 4.0 (compared to 2.5 obtained with the shear wave tomography approach). Nevertheless, the range in lateral viscosities here is rather narrow, in comparison to the lateral viscosity contrasts up to 10 4 found across southern British Columbia at 200 km depth (Figure 3 in Yousefi et al. (2021)).

Role of Material Compressibility
Previous models assumed incompressibility as well as laterally homogeneous viscosity. Compressibility leads to greater amplitude of predicted uplift rates in Southeast Alaska (Tanaka et al., 2015). However, differences between compressible and incompressible models can be reduced (to <10%) by adjusting the flexural rigidity of the elastic lithosphere (Tanaka et al., 2011). To isolate the role of compressibility we analyze the outcome of the best fit compressible 3-D model with its incompressible version, where the latter is obtained by adjusting Poisson's ratio to 0.4999. Note that this means that we only consider material compressibility (Klemann et al., 2003) and ignore the effects of internal buoyancy forces resulting from changes in material volume through deformation. The differences between the incompressible and compressible models can be seen in Figures 7c and 7f. The largest absolute differences, up to 5 mm/yr (∼15%), are seen in GB and YK ice fields, where the largest uplift signal is seen. However, the largest percentage change in rates (20%-25%) is seen in the Southeast corner of Figures 7c and 7f, corresponding to a lower uplift signal. The lithosphere in the compressible model can deform more easily than in the incompressible model, leading to higher peak uplift predictions. Figure 6a includes the misfit values of a select number of incompressible models. The best fitting incompressible model has a misfit 11.7% larger than the best compressible model, so the inclusion of compressibility is an important model element. The incompressible models require a lower background viscosity than the compressible models. In other words, the compressibility weakens the material, and we need to strengthen it to achieve the same uplift rate as the incompressible model by increasing the viscosity. The best fit averaged profile (the red line in Figure 5a) shows somewhat higher viscosity values (∼4.2 × 10 19 Pa-s) than the best fit model by  (3.0 × 10 19 Pa-s) in the asthenosphere.

Implications for Earth Structure
Upper mantle temperatures are widely considered responsible for most seismic velocity variations in the upper mantle, whereas compositional effects are thought to have second-order effects (e.g., Cammarano et al., 2003;. However, Goes and van der Lee (2002) also point out that the Western U.S. shows very low seismic velocities, which are likely due to fluids in the mantle introduced by the long history of subduction. In addition, Trampert and van der Hilst (2005) argue that chemical heterogeneity can introduce first-order uncertainties in the conversion from shear wave anomalies to temperature.
By using GIA and seismic tomography we were able to constrain the thermal effect on seismic velocity variations. Wang et al. (2008) showed that the thermal effects in Laurentia and Fennoscandia due to seismic anomalies are between 20% and 40%, by assuming a constant scaling factor throughout the mantle and not taking anelasticity effects into account. Wu et al. (2013) found a thermal contribution of 65% in the upper mantle, where they applied different scaling factors for the upper and lower mantle and anelasticity was taken into account.
We have found a low thermal contribution of 10% to the seismic anomalies, which implies that non-thermal effects control variations in seismic velocities. Although model misfit is similar if both background viscosity and thermal contribution are higher, such models lead to worsened fit relative to our preferred model (Figure 6a). GIA observables in Southeast Alaska are insensitive to the lower mantle structure and only weakly sensitive to deeper layers of the upper mantle . Therefore, it is not likely that scaling factors for deeper layers will have a large impact on our results. Uncertainty in  (Karato, 2008). Thus, the scaling factor may be underestimated between 10% and 20%. However, these uncertainties cannot explain the difference found with the global studies that found larger scaling factors.
Thus, we are left with the conclusion that non-thermal contributions to seismic velocities are large, and they are likely due to the presence of hydration and/or partial melt. If the mantle is substantially hydrated, anelastic effects would be stronger , and Equation 1 becomes larger, which results in a larger scaling factor. Considering that the region had a long history of past subduction, it is indeed likely that the mantle is substantially hydrated .

Model Limitations
There are a number of limitations to the GIA model that are briefly discussed. First of all, the uncertainties regarding the ice load model, both spatially and in time, influence the obtained earth model parameters. These uncertainties are related to both historic and PDIM load changes. The PDIM rates  were constrained by means of comparing GIA predictions with GPS observations in . Scaling of the ice thinning rates (Section 2.2) may not be uniform within Southeast Alaska and select areas may have an asynchronous ice load history with respect to the regional ice load model (such as YK and GB). The ice loading history was optimized for Southeast Alaska and this may not hold for all of Alaska. Moreover, the spatial loading history by Berthier et al. (2010) may be subjected to uncertainties and biases, as discussed in more detail in that study.
A second uncertainty relates to limitations in the earth model. The seismic velocity model used here (Schaeffer & Lebedev, 2013) has limited spatial resolution in this area because of the sparse station distribution. The assumptions in the seismic tomography inversion, such as model regularization, result in a very smooth model in this situation, so the details of the velocity structure are not well constrained. Regional seismic models exist (e.g., Jiang et al., 2018), but cover too small an area to have been used in this study. The lithospheric thickness was fixed for our main approach. It is possible that a good fit could obtained for a different combination of earth model parameters, for example a larger lithosphere thickness, larger background viscosity and larger scaling factor. We investigated an alternative approach based on an olivine flow law, in which lithospheric thickness is not specified a priori and the 3-D variations are derived from a global temperature model based on seismic and gravity data. However, it resulted in worse fits compared to the scaling of seismic anomalies because it effectively imposed a lithosphere that was much too thick. This is likely due to limitations in the data used in the global WINTERC-G model. Another limitation of this alternative approach is that only diffusion creep was modeled. A number of GIA studies have shown that a power-law rheology or composite rheology improved the overall fit to GIA observables (e.g., van der Wal et al., 2010). Viscosity required for best fit cannot be reconciled with the presence of hydration or standard grain sizes. Furthermore, background stresses could be significant here and other creep mechanisms, such as grain boundary sliding and transient creep could play a role. The latter could play an important role considering the timescales of the ice history. Transient creep has been shown to play a significant role in post-seismic studies in Alaska on monthly to decadal timescales . However, it is unknown how transient creep will affect the modeled uplift rates given the past and current ice load changes.

Conclusions
In this study, the shallow upper mantle viscosity structure beneath Southeast Alaska is studied using shear wave tomography and mineral physics in a GIA model for LIA and present-day ice thickness variations. The model is constrained by GPS uplift rates. The role of thermal effects on shear wave velocity anomalies is investigated by using an adjustable scaling factor, which determines what fraction of the seismic velocity variations is due to temperature changes, as opposed to non-thermal causes. If the scaling factor is 0, then there is no thermal contribution to variations in the seismic velocity. Contrarily, a scaling factor of 1 indicates that variations in the seismic velocity are only due to thermal effects. The viscosity values are computed using a law that scales seismic velocity anomalies, and relies on a background viscosity. The scaling factor also results in lateral variations in viscosity. By using this aspect, the role of lateral viscosity variations that are expected in this tectonically active region is also investigated. Our best fit model is obtained with a temperature scaling factor of 0.1 and background viscosity of 5.0 × 10 19 Pa-s. Models with a higher background viscosity and a higher scaling factor gave similar, but worse, misfit. An acceptable range of model parameters is 0.0-0.2 for the scaling factor and 3.7-7.0 × 10 19 Pa-s, which is based on an increase of misfit within 20% which approximates a 2σ uncertainty. This result implies that the contribution of thermal effects on shear wave velocity variations is small, which implies that seismic anomalies in the shallow upper mantle are mainly controlled by non-thermal effects such as hydration and/ or partial melt. The presence of hydration and/or partial melt  is consistent with the tectonic history of the region. For the best fit model, the viscosities at a depth of 80 km vary between 1.9 × 10 19 and 4.5 × 10 19 Pa-s and viscosity variations decrease within deeper layers. The viscosities obtained here are in the same range determined by previous 1-D GIA studies focused on Southeast Alaska Sato et al., 2011).
To address the relevance of lateral variations in the viscosity, we have compared the 3-D results to a radially symmetric model. The 1-D viscosity profile is obtained by averaging the 3-D viscosities in each Earth layer within a predefined area. The outcome shows that the 1-D model has a slightly better fit, however, the residuals cannot be distinguished from each other within measurement uncertainties of 2σ. Therefore, we conclude that 3-D variations do not have significant impact on the predicted uplift, given the current accuracy and spatial distribution of measurements.
Our best-fit viscosities are affected by uncertainty in the ice model and earth model as discussed in Section 3.7. Knowledge of the post-LIA ice history in the region is likely better than ice history in areas subjected to GIA due to post-LGM ice melting, and large improvements from field data are not expected. However, more detailed maps of current ice changes would lead to more accurate prediction of elastic component of the uplift rate. The conclusion that 3-D viscosity does not lead to a significantly better fit is mainly dependent on the accuracy of the seismic model. Although we have used the most recent data that is suitable for our model of the study region, more seismic data or improvements in the seismic inversion would change the pattern and amplitude of the anomalies. Furthermore, the conversion of seismic velocities to viscosity introduces further uncertainty. Apart from factors such as hydration and anelasticity, variation in the parameters used to compute the viscosity anomalies in Equation 1 are possible, see for example, Ivins et al. (2021) where uncertainty from these sources for Antarctica is quantified to be up to a factor of ∼3 (0.5 on a log-scale). Further progress in 3-D Earth rheology requires improvements in seismic models and laboratory deformation experiments to provide more accurate rheologic parametrizations that can be tested in geodynamic models.

Data Availability Statement
Input files of our FE models are available at 4TU Center for Research Data (https://data.4tu.nl/) under the name of this paper (https://doi.org/10.4121/16691329). GPS data are provided in . The shear wave velocity data we use are from the shear wave tomography model SL2013sv (Schaeffer & Lebedev, 2013). Most of the figures were prepared using Generic Mapping Tools (https://www. generic-mapping-tools.org/). review and constructive suggestions which helped us to improve this manuscript. Part of this work has been done in the framework of the project "3D Earth-A Dynamic Living Planet" funded by the ESA as a Support to Science Element (STSE). This research has been partly financially supported by the GO program of the Netherlands Organization for Scientific Research (NWO), project number: ALWGO.2019.001. J. T. Freymueller was supported by the NASA, award 80NSSC17K0566.