Long‐Wavelength Gravity Field Constraint on the Lower Mantle Viscosity in North America

The long‐wavelength negative gravity anomaly over Hudson Bay coincides with the area depressed by the Laurentide Ice Sheet during the Last Glacial Maximum, suggesting that it is, at least partly, caused by glacial isostatic adjustment (GIA). Additional contributions to the static gravity field stem from surface dynamic topography, core‐mantle boundary (CMB) topography, and density anomalies in the subsurface. Previous estimates of the contribution of GIA to the gravity anomaly range from 25% to more than 80%. However, these estimates did not include uncertainties in all components that contribute to the gravity field. In this study, we develop a forward model for the gravity anomaly based on density models and dynamic models, investigating uncertainty in all components. We derive lithospheric densities from equilibrium constraints but extend the concept of lithospheric isostasy to a force balance that includes the dynamic models. The largest uncertainty in the predicted gravity anomaly is due to the lower mantle viscosity, uncertainties in the ice history, the crustal model, the lithosphere‐asthenosphere boundary, and the conversion from seismic velocities to density have a smaller effect. A preference for lower mantle viscosities >1022 Pa s is found, in which case at least 60% of the observed long‐wavelength gravity anomaly can be attributed to GIA. This lower bound on the lower mantle viscosity has implications for inferences based on models for mantle convection and GIA.


Introduction
The global gravity model XGM2016 exhibits a negative anomaly of about 50 mGal near Hudson Bay for wavelengths larger than 600 km ( Figure 1) (Pail et al., 2018). The shape of this anomaly resembles the contours of the deflection due to the former Laurentide Ice Sheet. Hence, the anomaly is thought to be caused by the incomplete rebound following the deglaciation of the Laurentide Ice Sheet (Kaula, 1972;Walcott, 1973), a process known as glacial isostatic adjustment (GIA). Because of incomplete GIA, the topography is not in equilibrium, and this topographic deflection can be seen in the gravitational potential field, hereafter referred to as the gravity field, or rather in anomalies of this field with respect to the gravity field of a reference Earth for which we use gravity anomalies. If GIA were the only cause for the gravity anomaly, the observed gravity anomaly with its small error would form a useful constraint on GIA models. However, the gravity field contains contributions from the crust, the lithosphere, the mantle, and the GIA. Before using the gravity anomaly to constrain GIA, these contributions need to be quantified.
The crustal and lithospheric contributions to the gravity field are from density contrasts at geological boundaries. One of these density contrasts marks the boundary between the crust and the mantle, the Moho. Knowledge of the geometry of this boundary is therefore important for gravity modeling. A second boundary is the lithosphere-asthenosphere boundary (LAB), which can be defined as a boundary separating the conductive and convective regimes (e.g., Eaton et al., 2009;Fischer et al., 2010;Sleep, 2005). This boundary is not characterized by a large jump in density but determines where the mantle can start to flow to equalize the weight of the overlying material. The LAB is therefore an important boundary for the force balance in the upper mantle and can be inferred from estimates of, among others, heat flow or seismic tomography Eaton et al., 2009). Beneath Hudson Bay, the lithosphere is cratonic and has a thickness of 150-200 km (Eaton & Darbyshire, 2010).
Another factor determining the gravity field contribution of the crust and the lithosphere is isostasy. Isostasy implies that the pressures at a certain depth are equal. For example, the crustal thickness, which delivers the buoyancy to maintain the topography, can be determined in the classical Airy isostasy theory (crustal isostasy Watts, 2001). In other studies, isostasy is calculated based on a lithosphere floating on top of a homogeneous asthenosphere (lithospheric isostasy Lachenbruch & Morgan, 1990). Isostasy is often a necessary assumption to be able to fit the observed gravity anomaly, also for North American studies (e.g., Métivier et al., 2016). However, lithospheric isostasy neglects forces from the mantle. Here, we will extend the lithospheric isostasy to include forces from GIA and convection models, which are expected to play a significant role in North America.
Mantle contributions to the gravity field consist of (i) density anomalies in the mantle, (ii) surface dynamic topography (Hager et al., 1985), and (iii) topography of the core-mantle boundary (CMB). (ii) and (iii) depend on the viscosity contrast between upper and lower mantle; a smaller contrast results in a larger signal for the surface dynamic topography, while the signal is smaller for the CMB topography. (i) and (ii) have opposite signs; a positive density anomaly drags the surface down, resulting in negative dynamic topography, which counteracts the positive gravity anomaly from the density anomalies. The lower boundary of the mantle marks the largest density contrast in the Earth, larger than the density contrast at the surface. Only long-wavelength surface features are affected by CMB topography, and for these wavelengths, it is hence important to include (iii). The long-wavelength signal in the gravity field and the geoid can be matched well by mantle convection modeling using seismic tomography as input for the mantle density distribution (Hager et al., 1985). For North America, the main mantle signal that is expected is that of the subducted Farallon slab, of which the geometry and subduction history are complex (Sigloch, 2011).  (Pail et al., 2018) over North America up to degree 60 (a) and up to degree 15 (b). The study area is indicated by the dashed red line, which contains the points where the value is at least 50% of the peak value.
The GIA contribution to the gravity field is mostly dependent on the ice sheet history and on the viscosity of the Earth's mantle. The ice thickness controls the magnitude of the deflection, while the viscosity controls how fast the equilibrium is reached. A large viscosity leads to a smaller displacement during the loading phase, as compared to low viscosities, and also to a slower rebound after removal of the load. Meanwhile, in a low viscosity mantle, the deflection during loading is larger, but a fast relaxation after unloading can also lead to a small remaining displacement. Because of this, there are often multiple solutions to the viscosity for the same displacement. This explains part of the differences in viscosity profiles obtained in different studies.
Previous studies attribute different percentages of the gravity anomaly to GIA (e.g., Métivier et al., 2016;Paulson et al., 2007). This discrepancy can, to a large extent, be explained by different assumptions of the underlying mantle viscosity and whether GIA, surface dynamic topography and crustal and mantle density variations are considered in the modeling or data correction. The first studies that try to explain the gravity anomaly over Hudson Bay note that it cannot be explained by GIA using a lower mantle viscosity of 10 21 Pa s and suggest by inference that the major contribution is that of mantle convection (Cathles, 1975;James, 1992;Peltier et al., 1992). Simons and Hager (1997) find that the GIA contribution is significant and that about 50% of the free-air gravity anomaly can be explained by GIA. In their study, they employ a lower mantle viscosity that is close to 10 22 Pa s. Tamisiea et al. (2007) used time-variable gravity from the GRACE mission to isolate the GIA signal and found viscosities between 10 21 and 10 22 Pa s. Consequently, they attribute only 25-45% of the free-air gravity anomaly to GIA. Finally, Métivier et al. (2016) investigated gravity, together with gravity gradients, by combining a lithospheric model with models for GIA and mantle convection and found values of at least 10 22 Pa s. In their study, GIA contributes more than 80%. All in all, the contribution of GIA to the gravity anomaly is still uncertain, with most recent estimates ranging from 25-45% (Tamisiea et al., 2007) to more than 80% (Métivier et al., 2016), with part of the spread explained by the unknown mantle viscosity.
The mantle viscosity is not well constrained, and many studies have attempted to determine its value by employing constraints on mantle convection models (e.g., Soldati et al., 2009;Steinberger & Calderwood, 2006), GIA models (e.g., Paulson et al., 2007;Wu & Peltier, 1983), or a combination of both (Forte & Mitrovica, 1996;Mitrovica & Forte, 2004). Mantle convection studies that determine the viscosity are often global studies employing a radial viscosity profile containing many layers. These global inferences of the viscosity generally favor lower mantle viscosities >10 22 Pa s (e.g., Bower et al., 2013;Perry et al., 2003;Steinberger & Calderwood, 2006). GIA studies have been performed on both a global and regional scale and generally use relative sea level (RSL) data and geodetic data to constrain the viscosity. No consensus has been reached on the value of the viscosity in the lower mantle under North America, with values varying by 2 orders of magnitude (10 21 to 10 23 Pa s) (e.g., Métivier et al., 2016;Paulson et al., 2007).
Dynamic models (i.e., the mantle convection model and the GIA model) contain more uncertain parameters than the viscosity. One of these parameters is the ice history, to which the gravity field is especially sensitive at the margin of ice sheets (Mitrovica et al., 1994;Wu, 2006). However, the extent of the Laurentide Ice Sheet is relatively well known, in contrast to its thickness (Stokes, 2017). Uncertainty due to the ice history has been included in some previous studies of the gravity field (James, 1992;Métivier et al., 2016), but not in all (Peltier et al., 1992;Tamisiea et al., 2007). For mantle convection modeling, density anomalies are needed, which are commonly derived from seismic velocity anomalies. The conversion factor between velocity anomalies and density anomalies can vary between 0.2 and 0.4 (Karato, 1993;Trampert et al., 2004), as determined by measurements and employed in convection models (Steinberger & Calderwood, 2006). This conversion factor can have a large influence on the resulting gravity, as it can amplify or minimize gravitational signals from the mantle. Métivier et al. (2016) assign a conversion factor to each viscosity layer in their mantle models but do not show the sensitivity to this parameter.
We combine contributions from a crust-lithosphere-asthenosphere (CLA) model, mantle densities, and CMB topography in a forward model and vary parameters of all components. If all geometries and densities were perfectly known, this would be able to explain the observed gravity anomaly. However, the assumption of a force balance is required to compute the densities in the CLA model. We include GIA and surface dynamic topography in the force balance for North America, which we argue is the way to include all model components consistently.
Viscosity controls the contribution from the dynamic models to the force balance and the contribution of CMB topography. In the presence of other uncertainties, viscosity is the most important parameter, and we are able to constrain it from fitting the observed gravity field. From this viscosity, we obtain the contributions of GIA and mantle convection to the gravity anomaly.
Section 2 explains the approach used to construct the density model of the CLA and its conversion to gravity anomalies. After that, we elaborate on the GIA and mantle convection models used and how they are included in the isostatic balance of the CLA. Section 3 starts by investigating the effect of uncertainty of the crustal and lithospheric model on the modeled gravity field. Next, we show the GIA and mantle contributions to the gravity field as a function of the viscosity profile. After this, we discuss uncertainties due to the ice history and due to the conversion from seismic velocities to density. We find the best-fitting solution for an earth model with varying upper and lower mantle viscosities and obtain a lower bound on the lower mantle viscosity.

Methodology
In this section, we start by stating the complete model used in this study. The top layers are the crust, lithosphere, and asthenosphere. Second, we explain how GIA and surface dynamic topography are incorporated in our CLA model using isostatic equilibrium. After that, we elaborate on the CLA models, the mantle models, and the GIA models. It is not our aim to explain small-scale features in the data, and we would not be able to do so with the current models. We therefore set a maximum spherical harmonic (SH) degree for the data-model comparison and explain the reason for this choice of maximum degree at the end of section 2.6.

The Complete Model
Our forward model includes the gravity field (Δg) of a CLA model, mantle density anomalies below 300 km (ρ m ), and topography at the CMB: Note that the signal of GIA and surface dynamic topography do not appear explicitly in Equation 1. These nonisostatic pressures are included in the isostatic balance of the CLA model (section 2.2). The mantle density anomalies are derived from seismic velocity anomalies; see section 2.4. The gravity signal from the topography at the CMB is computed by a mantle convection model (Tosi, 2008).
It might seem contradictory that the influence of GIA is not explicit in Equation 1, yet we are able to constrain its contribution to the gravity field. The explanation is that the GIA contribution is determined by viscosity, which controls the contribution of the mantle in the last term of Equation 1 but also the contributions of GIA and surface dynamic topography in the force balance that is employed for the CLA model as will be explained in section 2.2. Thus, a data-model comparison allows us to constrain the viscosity, and with the viscosity, we are able to retrieve the contribution of GIA and mantle convection to the gravity field.
For Equation 1, the gravity anomaly needs to be computed for layers of arbitrary density. To do this, we use a spectral method that transforms 3-D spherical density models into SH coefficients of the gravitation potential (Root et al., 2016) and employ the SHTools package (Wieczorek & Meschede, 2018) to convert the coefficients to gravitational potential fields in the spatial domain. We calculate the gravity anomaly at the height of 6,738 km, the height at which XGM 2016 is calculated.
Although the geoid is commonly used, we opt to show gravity anomalies, which is the radial derivative of the gravity potential. Our choice to represent the gravity field is, in principle, arbitrary for the long wavelengths that we employ in this study. Strictly speaking, we are computing gravity disturbances. These disturbances are referred to as gravity anomalies in this study.

Isostasy in the CLA Model
In principle, the gravity field can be represented by geometry and density information of each layer in the subsurface. In practice, accurate density information is not available at each depth, and the assumption of isostasy is made to solve for densities or geometry. In this study, we employ lithospheric isostasy (Lachenbruch& Morgan, 1990;Root et al., 2017), which involves adjusting the density of the lithosphere. In an isostatic balance, no dynamic processes should be present. However, the dynamic forces of GIA and surface dynamic topography are acting on the bottom of our CLA model, so we need to correct for those. We present a consistent way of including information from crustal models, isostasy, and GIA and surface dynamic topography models.
To implement lithospheric isostasy, free body diagrams are made for mass columns up to 300 km ( Figure 2). The forces involved are those caused by the weight of the crust, lithospheric mantle, and asthenosphere, and we implement GIA and surface dynamic topography as radial forces acting at 300 km. The pressure at 300-km depth for each column should equal that exerted by a reference column. Here, the reference column consists of a 30-km-thick crustal layer and a 270-km-thick mantle layer with densities of 2,850 and 3,300 kg/m 3 , respectively ( Figure 2). Equilibrium of the forces is then achieved in the following manner: For each layer the pressure at 300 km can simply be calculated from its weight per area. To include surface dynamic topography, we calculate the radial stress caused by this process (Flament et al., 2013), and add it to the force balance of Equation 2. The density (ρ asth ) is that of the asthenosphere and equal to 3,300 kg/m 3 , g is the average value of the gravity, and h DT is the height of the surface dynamic topography. σ rr is calculated by the mantle convection model (Tosi, 2008), discussed in section 2.4.
In principle, the contribution of GIA above 300 km is already included in the geometry of the crust and lithosphere because their boundaries are deflected by GIA. However, GIA below 300 km plays a role through the stress that acts at the bottom of our CLA model, similar for the surface dynamic topography. These stresses affect the isostatic balance, so we include stresses from both GIA and surface dynamic topography in the isostatic balance to create a force balance (Equation 2). Following the approach of Root et al. (2015), we implement the effect of GIA by shifting the layers above 300 km according to the respective GIA deflection at the surface, h GIA , defined as positive upward. This way, we assume isostasy based on a configuration in which GIA is no longer present. h GIA is calculated by the GIA model, discussed in section 2.5, and is taken to be the same for all elastic layers, which are the lithospheric layers that conform to Root et al. (2015).
Combining these ideas in the force balance (Equation 2), assuming constant gravity in the top 300 km yields The summation i is over all the crustal layers defined in the crustal model, where r ub,i and r lb,i are the upper and lower boundary of layer i, respectively, and ρ crust,i is the crustal density for the layer, which can vary laterally. An earlier version of this equation, without the processes of GIA and surface dynamic topography, is The forces exerted by the crustal, lithospheric, and asthenospheric layers are denoted by F crust , F litho , and F asth , respectively. F GIA and F DT are the forces due to GIA and surface dynamic topography.
shown in Root et al. (2017). The first, second, and third terms on the left-hand side of the equation represent the crustal, lithospheric, and asthenospheric layers in the model, respectively, and the right-hand side represents the reference column. r 0 is equal to 6,371 km, the radius of the Earth. The radii to the Moho and the topographic boundaries are defined as positive upward. We assume that the geometry and density of the crust are reasonably well known from seismic data compared to deeper layers. However, deriving density from seismic models in the asthenosphere does not result in a good fit with gravity data. Therefore, isostasy is used as a constraint to achieve a density model that provides a good fit with gravity. From this constraint, we can only constrain one parameter. Even if there were reliable density information from seismic models, there would be a strong trade-off effect between different LAB or asthenosphere density and obtained lithospheric density. In the end, we opt to adjust the density of the mantle lithosphere, which is less well known and is represented by Δρ in Equation 4. It is important to note that the GIA contribution to the first four terms contains the entire GIA contribution (Root et al., 2015). Boundaries below 300 km have a smaller density change and/or a smaller deflection, and these are neglected. Similarly, the fourth term contains all of the effects of surface dynamic topography. This term is negative because σ rr is defined as positive upward in Equation 3, and, consequently, the direction of this load is opposite that of the gravitational loads. Thus, a positive surface dynamic topography contribution results in an effective negative mass that will be compensated because of the pressure balance represented by Equation 4. To recapitulate, the force balance (Equation 2) is transformed into a pressure balance. We can do so because the area of the model columns is the same as that of the reference column. Assuming that gravity is constant for all layers in the CLA model, we obtain the pressure balance shown in Equation 4.

Crust, Lithosphere, and Asthenosphere Models
The CLA model relies on seismic information for crustal thickness, crustal densities, and a LAB estimate, together with the lithospheric densities obtained in section 2.2. To account for uncertainty in the crustal thickness, two crustal models are used: CRUST1.0 (Laske et al., 2013) and a crustal model based on the U. S. Geological Survey (USGS) Global Seismic Catalog (GSC) database, which was interpolated to a 1 × 1°g rid using kriging interpolation . This data set has been augmented over North America with data from the Geological Survey of Canada (Schetselaar & Snyder, 2017) and will be named GSCaug hereafter. CRUST1.0 has a resolution of 1×1°, and each cell has a unique eight-layer density profile. The GSCaug data set only presents the Moho depth. We adopt a crustal density of 2,850 kg/m 3 for the GSCaug data set, based on the reference profile described in section 2.2. For both crustal models, the topography, bathymetry, and ice-cover are taken from CRUST1.0, as uncertainties in these components are negligible for the long-wavelength signal studied in this article. Moho depths of the crustal models are shown in Figures 3a and 3b. In oceanic areas, the Moho depth is 20 km at most. The Moho depth is clearly larger for continental areas, with values of 30 to 50 km. Around Hudson Bay, there are regional differences of up to 10 km between the crustal models. CRUST1.0 is used as the default crustal model for the rest of the analysis. This is different from Métivier et al. (2016), in which an isopycnal configuration of the crust was used.
The lithosphere and the asthenosphere are separated by the LAB. To account for uncertainty in this depth, we use two estimates for the LAB (Figure 3c and 3d). The first option is the LAB model of Hamza and Vieira (2012), a model that is derived from estimates of surface heat fluxes and crustal modeling. The second option is obtained from WINTERC v5.2, which is a 3-D model of the lithosphere and the upper mantle based on a joint inversion of surface waveform tomography, surface heat flow, and elevation of the topography. In both models, the LAB reaches its largest depth in an area below Hudson Bay, thereby correlating with the observed gravity anomaly. In the LAB model compiled by Hamza and Vieira (2012), the LAB low over Hudson Bay is more confined and larger in amplitude than the WINTERC LAB.
We have used the LAB from Hamza and Vieira (2012) as our reference model and the LAB estimate from WINTERC v5.2 as our alternative model in the remainder of this study. We have to note that the LAB from Hamza and Vieira (2012) is probably not well constrained since estimates of surface heat fluxes form a poor constraint in terms of sparsity and error. However, since there is no real density jump at the LAB, we do not expect large changes in the gravity field as a result of choice here. The sensitivity to this choice is investigated more thoroughly in section 3.

Mantle Below 300 km
The contributions of the mantle below 300 km are computed using a mantle convection model (Tosi, 2008). This spectral finite element code solves the incompressible Stokes problem and computes the gravity field resulting from density anomalies and boundary deflections. For the radial direction, the model employs finite elements, while for the angular direction, SHs are used to parameterize the solutions to the Stokes problem. The model uses mantle density anomalies as input and produces dynamic topography at the surface (Hager et al., 1985) and at the CMB. It is assumed that convection is only driven by density perturbations below 300 km.
The density anomalies are, in turn, derived from seismic velocity anomalies. Here, the seismic velocity anomalies are taken from the global, composite tomography model SMEAN2, which uses the approach of Becker and Boschi (2002) to combine S40RTS (Ritsema et al., 2011), GyPSUM-S (Simmons et al., 2010), and SAVANI (Auer et al., 2014). There is uncertainty inherent in the choice of tomography model. However, SMEAN2 is an aggregate of multiple tomography models, and Becker and Boschi (2002) show that the tomography models correlate for long wavelengths. Moreover, SMEAN2 is an aggregate of multiple tomography models. We assume the uncertainty is smaller than the uncertainty in conversion factor (Equation 5) and that the effect is partly compensated by varying the conversion factor. Shear-wave velocity anomalies (Δv) can be converted to density anomalies (Δρ) by a conversion factor (p) (Karato, 2008): In this study, the conversion factor has a constant value of 0.15. In reality, the conversion factor can change radially (Karato, 2008;Steinberger & Calderwood, 2006) and laterally, but a single value is sufficient if the sensitivity to the parameter is small. The uncertainty introduced by the conversion factor is analyzed in section 3. We do not include a transition zone due to trade-off effects, which allow both low and high viscosity in the transition zone to fit the geoid (King, 1995). The layering is the same as the GIA model, in which trade-off effects are perhaps even stronger (Paulson et al., 2007). The converted density anomalies, together with a three-layered viscosity profile (elastic lithosphere, upper mantle, and lower mantle), are used as input in the mantle convection code. The lithosphere is assumed to have a thickness of 100 km. In Appendix A1, we investigate the sensitivity to the lithospheric thickness in the mantle convection model. The viscosities of the upper and lower mantle are separated by the 670 km discontinuity. The density of the core is assumed to be homogeneous and is set equal to 4,500 kg/m 3 . The CMB and the Earth's surface are modeled as free-slip, impermeable boundaries. The output of the mantle convection code consists of stresses at the top and bottom boundaries. The stresses are converted to surface dynamic topography values with Equation 3. After, the surface dynamic topography is converted to pressure to compute lithospheric density anomalies in order to fulfill isostatic balance in the CLA model using Equation 4. The mantle density anomalies and the CMB topography are converted to SH coefficients following the approach of Root et al. (2017), and these coefficients are added to the coefficients from the CLA model, as shown in Equation 1.

GIA Models
The GIA response to the glacial loading is calculated with a normal mode method (Wu & Peltier, 1982) for radial variations in viscosity within a multilayer model (Vermeersen & Sabadini, 1997) with self-consistent sea levels (Mitrovica & Peltier, 1991c). Rotational feedback (Milne & Mitrovica, 1998;Wu & Peltier, 1984) and geocenter motion (Greff-Lefftz & Legros, 1997) are incorporated in the model. Since we will only look at degrees 2 to 15, as will be shown in section 2.6, the direct effect of geocenter motion will drop out, but any coupling to higher terms through the sea-level equation, although small, is included. The code is developed by Schotman (2008) and has recently been benchmarked for simple loading scenarios in Martinec et al. (2018). In general, viscosity in GIA models can vary radially (1-D) or in both the radial and lateral directions (3-D). GIA models with a 1-D viscosity in North America match results from laterally averaged 3-D models reasonably well (Geruo et al., 2013). Moreover, the effect of 3-D viscosity on predictions around Hudson Bay is limited (Li et al., 2020). Therefore, it is expected that the 1-D Earth model employed in our model produces reasonably accurate results.
The GIA model adopts a similar three-layer Earth model as the mantle convection code discussed in section 2.4, consisting of an 80-km-thick, elastic lithosphere and a viscous upper (<670 km) and lower (>670 km) mantle. Two viscous layers in the mantle provide sufficient degrees of freedom to model GIA observables (Paulson et al., 2007). Therefore, we do not consider separate layers to model shallower phase transitions. The elastic parameters are obtained from the preliminary reference Earth model (PREM Dziewonski & Anderson, 1981) and are the same as in van der Wal et al. (2009). This lithospheric thickness is different from that of the mantle convection model, but since our results turned out to be insensitive to the lithospheric thickness, this difference will not have a large effect.
An important uncertainty in the GIA model is caused by the unknown ice loading history. Four different ice histories are employed to assess this uncertainty. The ice models are ICE-6G (Argus et al., 2014;Peltier et al., 2015), the model by Lambeck et al. (2017), which will be labeled LW-6, and two variants of the GLAC-1D model (Tarasov et al., 2012), named GLAC-1D nn9894 and GLAC-1D nn9927. The ICE-6G and GLAC-1D models are global models, while LW-6 is a regional model. ICE-6G uses ice extent constraints and is tuned to fit RSL data and geodetic constraints, although the fitting started with a model based on ice dynamics. The North American sector of GLAC1-D uses much of the same RSL and geodetic constraints as that of ICE-6G, as well as marine limit and strandline data. It also accounts for age uncertainty in the geologically inferred deglacial margins and is derived from an approximate Bayesian formalism applied to a thermo-mechanically coupled glaciological model. The ICE-6G and GLAC-1D models are based on the VM5a viscosity profile (Peltier et al., 2015), while the viscosity profile used for the LW-6 model consists of a three-layer viscosity profile with an upper mantle viscosity of 5.1 × 10 20 Pa s and a lower mantle viscosity of 1.3 × 10 22 Pa s. There will be a bias in the viscosity fit toward the implicit viscosity profiles, although gravity data have not been used to constrain the ice models. To reduce the effect of the bias, we use the different ice models and especially the GLAC-1D model, which depend less strongly on the viscosity because they are primarily controlled by ice dynamics.
Ice thicknesses at 26 ka are up to 5,000 m in the ICE-6G model and up to 4,000 m in the other models ( Figure 4). ICE-6G also has thicker ice in the western part of North America compared to the other models. Together, this is a partial representation of the uncertainty in the ice loading history. For all models, three glacial cycles of 112 ka are used to account for the effect on the gravity anomaly of earlier glaciations in models containing larger values for the lower mantle viscosity. The glaciation phase, between 26 and 112 ka, of LW-6 is not available. We use the ICE-6G ice model to fill this gap. The first two glacial cycles are assumed to be the same as the last one.

SH Truncation Limit
The signal that we want to explain is the long-wavelength gravity field, which contains most of the GIA and mantle convection. The truncation limit should be a trade-off between containing most of the GIA and mantle signal and minimizing the uncertainties in the other components, especially in the crustal model, which can introduce uncertainties up to 110 mGal (Root et al., 2015). Also, the lithospheric density anomalies are especially uncertain in the short-wavelength region. Another argument in favor of a low maximum SH degree is the assumption of local isostasy made in the model, which works best for long-wavelength signals (Gvirtzman et al., 2016;Watts, 2001), since flexural isostasy starts to contribute significantly to degrees larger than ∼30 (Watts & Moore, 2017).
Mantle convection manifests itself in longer wavelengths and contains most of its signal below SH degree 10 (e.g., Gu et al., 2001;Steinberger et al., 2019;Su & Dziewonski, 1991). Therefore, the truncation is mostly determined by the GIA signal. In Figure 5, the amplitude and the location of the GIA signal are plotted for models containing an upper mantle viscosity of either 2 × 10 20 or 4 × 10 20 Pa s and a lower mantle viscosity >10 21 Pa s. The solid lines result from models with an upper mantle viscosity of 4 × 10 20 Pa s and a lower mantle viscosity of 3.2 × 10 21 Pa s (blue), 1.3 × 10 22 Pa s (green), and 2.6 × 10 22 Pa s (red). The gravity anomaly for different viscosity profiles is represented by the shaded areas and exhibits the same behavior as shown by the solid lines.
The idea is to find a truncation limit above which the GIA gravity signal loss is relatively small. The amplitude of the GIA signal over North America starts to decrease slightly for a maximum SH degree lower than 20, stabilizes again, and then decreases rapidly for a maximum SH degree lower than 10 (Figure 5a). A second criterion is based on the location of the maximum amplitude of the GIA signal in the models. For different truncation limits, the location of the maximum amplitude in the gravity field is compared to that of the original model, which uses a maximum SH degree of 50 (Figure 5b). For a truncation limit at SH degree 10, the distance to the original maximum amplitude is almost 300 km. The previous sentences argue for a cut-off that is higher than 10 but does not have to be higher than 20 to capture most of the GIA model. Another criterion is that we want to reduce the effect of uncertainty in the crustal model. Since the large-scale features of the crustal models are very similar for a maximum SH degree of 15 or less, and the uncertainty in the crustal model is increased significantly for SH degree 20 ( Figure C1), we will use SH degree 15 as the maximum degree for all models and observations in the rest of this study. We also briefly looked at cut-off degrees 14 and 16.

Results
To obtain density anomalies that are necessary to forward model gravity, we use the force balance of Equation 4. We fit the viscosity to observations, and from our preferred viscosity, we obtain the GIA and mantle convection contributions to the gravity field. As mentioned in section 2.2, the lithospheric mantle densities are adjusted to ensure a force balance. We assign a single lithospheric density to each grid cell between the Moho and the LAB. We obtain lithospheric densities between 3,320 and 3,380 kg/m 3 ( Figure B1a) using the LAB from Hamza and Vieira (2012). The range in lithospheric densities is larger (3,140-3,500 kg/m 3 ) for the WINTERC v5.2 LAB ( Figure B1b). The largest differences between the modeled densities are present in areas where the lithosphere is relatively thin since here the densities need to be adjusted more to accommodate a similar change in LAB. In areas where the lithosphere is thick, like Hudson Bay, differences in the modeled lithospheric densities are less prominent. Thus, because the LAB is used to determine the lithospheric density needed for isostasy, the sensitivity of the gravity anomaly to the LAB estimate is reduced. Figure 6a shows the gravitational signal due to a combination of our crustal model (CRUST1.0) and our LAB (taken from Hamza & Vieira, 2012). A small gravity low of up to 15 mGal is present just Southwest of Hudson Bay. This gravity anomaly extends to the south and reaches 30 mGal south of Lake Michigan. The gravity high over the Rocky Mountains is up to 20 mGal. The uncertainty due to the crust can be caused by (i) the density profile adopted and (ii) the Moho employed in our model. To determine the effect of different density profiles, we compare the gravitational signal from a layered density profile with that of an isopycnal crust with a density of 2,850 kg/m 3 without changing the Moho (in both cases, the Moho is that of 10.1029/2020JB020484

Journal of Geophysical Research: Solid Earth
CRUST1.0). The fact that we only focus on SH degrees 2 to 15 greatly reduces the range, as only the long-wavelength signals remain, and these are generally more consistent among different crustal models (Figures 6b-6d). The range in the crustal signal over Hudson Bay is small and, for the most part, below 5 mGal. In two regions, the range is 10 mGal, namely, in the Canadian Arctic Archipelago and in the geologically complex Rocky Mountains. Figure 6c shows the range in the crustal signal due to the Moho. This is the spread in gravity signal caused by employing the CRUST1.0 and the GSCaug Moho models. When determining the range caused by different Mohos, we have made use of an isopycnal crust. The range in the signal due to the Moho is small and, with the exception of the region over the Canadian Arctic Archipelago, below 5 mGal. We determine the range in signal due to different LABs in the same way as we did for the Moho. The range in the gravity contribution due to different LAB representations is only up to 5 mGal over Hudson Bay (Figure 6d). The reasons for these small numbers, despite large differences in LAB, are the compensating effect of fitting lithosphere densities to the isostasy constraint and the absence of a density jump at the LAB.
GIA and surface dynamic topography both contribute to the total modeled gravity signal of the CLA model through their contribution to the force balance (Equation 2) and the resulting effect on the lithospheric density. To compare the contributions of GIA and the mantle, the GIA and surface dynamic topography heights calculated in section 2 are converted to SH coefficients. The contributions from mantle density anomalies and CMB topography are added to that of surface dynamic topography to form the total effect of the mantle. From the SH coefficients for GIA and the mantle, the gravity anomalies can be calculated and compared. We vary the viscosity values of the upper and lower mantle between 10 20 and 10 23 Pa s and calculate the gravity signal at the location of the modeled minimum in the gravity anomaly. We do so because we believe that the gravity anomaly minimum characterizes the signal better than a value of the gravity anomaly at a fixed location. The viscosity influences the gravity signal through the GIA and surface dynamic topography contribution to the force balance and also through the CMB topography. Thus, even though the GIA deflections can not be isolated from the crustal geometry, we can obtain the GIA contribution to the gravity field with a

10.1029/2020JB020484
Journal of Geophysical Research: Solid Earth strong constraint on lower mantle viscosity, with which we can predict the GIA signal. Figure 7 exhibits a wide range of values, depending on the viscosity profile. The GIA contribution is most sensitive to the viscosity of the lower mantle. For lower mantle viscosities >10 22 Pa s, GIA contributes at least 30 mGal to the negative anomaly. This contribution decreases when the lower mantle viscosity decreases. Lower viscosities imply a shorter relaxation time, and consequently, less remaining uplift is present in the lithosphere. This results in a smaller contribution to the gravity field. For all viscosity profiles, the contribution of the mantle below 300 km (including surface dynamic topography) does not exceed −20 mGal and can even be weakly positive for lower mantle viscosities >10 22 Pa s. The crust and lithosphere contributions do not depend on the underlying viscosity profile and contribute 15 ± 10 mGal to the gravity anomaly ( Figure 6a). The spread in the crustal signal comes from the changing location of the minimum and the use of different crustal models.
As we have seen that the gravity field is sensitive to the lower mantle viscosities, it is insightful to exhibit the contributions of GIA and the mantle below 300 km (with surface dynamic topography) and their dependence on the lower mantle viscosity. Depending on the lower mantle viscosity, the GIA signal can be up to 20 mGal (Figure 8a) or up to 40 mGal (Figure 8b). For mantle convection, amplitudes are lower, but the sign can be reversed depending on the lower mantle viscosity. For a lower mantle viscosity of 3.2 × 10 21 Pa s, the signal is weakly negative and consequently contributes positively to the observed gravity anomaly (Figure 8c). For a viscosity of 2.6 × 10 22 Pa s, anomalies due to mantle convection are weakly positive over Hudson Bay (Figure 8d). The positive signal in Figure 8d compensates slightly for the increased amplitude of the negative anomaly due to GIA.
Next, we look for constraints on the lower mantle viscosity, taking into account uncertainties in other components. In order to place constraints on the viscosity of the lower mantle, we have created models according to Equation 1 for different combinations of the upper and lower mantle viscosity, using the crustal and lithospheric models of Figure 3. We compared the results of these models with the observed gravity field for the area depicted by the red dashed line in Figure 1, which is the area covering Hudson Bay and the region south of Hudson Bay up to major lakes like Lake Michigan. A misfit is then calculated using the following formula: where N is the number of grid cells inside the area of the latitude/longitude grid, and o i and p i are the observations and the predicted values at gridpoint i, respectively. σ i is the uncertainty of the observation.
Since the values of the observed gravity anomaly at each gridpoint are correlated, we cannot use this as the uncertainty. Instead, we attribute an arbitrary value of 1 to σ i and calculated χ 2 values for four possible combinations of crustal models and LABs. We then averaged the results of these models. The relative value with respect to other viscosity profiles is what matters in these plots. We have highlighted the

10.1029/2020JB020484
Journal of Geophysical Research: Solid Earth models for which χ 2 is below 40. This is an arbitrary boundary created to emphasize better-performing models. The misfit for different upper and lower mantle viscosities is shown in Figure 9. The well-performing models are found almost exclusively for lower mantle viscosities of more than 10 22 Pa s ( Figure 9). Models containing lower mantle viscosities in the range 10 21 to 10 22 Pa s underestimate the negative anomaly in the gravity field observed over Hudson Bay, naturally resulting in high χ 2 values.
Next, we investigate whether this conclusion changes for different choices in the modeling. The good fit for lower mantle viscosities above 10 22 Pa s does not change if we change the lithospheric thickness in the GIA model from 80 to 115 km or to 150 km ( Figure 10) or if we make the area of interest smaller (the area delimited by the dashed red line in Figure 1). For example, we can choose the area to contain all points that have a value that is at least 40% of the peak value, as opposed to the 50% threshold used in the rest of this study. Moreover, if simulations were performed with a different SH truncation limit (e.g., 14 or 16 as the upper limit), the general patterns in the misfit plot remain the same. The most important other sources of uncertainty are discussed in the following paragraphs.
The next parameter that we will discuss is the ice history, which is used as an input to the GIA model. Variations in ice heights and the time of melting translate directly in the gravity signal (Mitrovica & Peltier, 1991b). The subplots in Figure 9 correspond to the four ice histories used. Lower mantle viscosities >10 22 Pa s show lower misfit values, regardless of the ice model used. This confirms that the preferred viscosity profile does not depend strongly on the ice history. For the GLAC-1D nn9894 ice history, lower mantle viscosities of 6.4 × 10 21 Pa s also perform well. However, for these specific well-performing models, the upper mantle viscosity needs to be >10 21 Pa s, which is not corroborated by other studies on the viscosity of the upper mantle in North America (e.g., Paulson et al., 2007;Sasgen et al., 2012;Tamisiea et al., 2007;Wolf et al., 2006). Thus, regardless of the employed ice history, lower mantle viscosities >10 22 Pa s are preferred.
The final parameter that we will test the sensitivity to is the conversion factor from seismic velocity anomalies to density anomalies. We vary the conversion factor between 0.1 and 0.4 to represent the range of possible values (Trampert et al., 2004) and study its effect on our conclusions. For each conversion factor, a χ 2 misfit is calculated using Equation 6. The largest sensitivity is to the viscosity of the lower mantle. For this The viscosity of the upper mantle is the same across all subplots and equal to 4 × 10 20 Pa s.

10.1029/2020JB020484
Journal of Geophysical Research: Solid Earth reason, Figure 11 shows the spread in χ 2 values as a function of the lower mantle viscosity, while the upper mantle viscosity is kept fixed at 4 × 10 20 Pa s. Almost all models containing lower mantle viscosities >10 22 Pa s have a lower χ 2 value than models containing lower mantle viscosities <10 22 Pa s, independent of the conversion factor used. The spread in misfit values between observations and models decreases when the lower mantle viscosity is increased. This is because, for those viscosities, the contribution of the mantle convection signal is close to zero or just about positive over North America (see Figure 7). Consequently, an amplification or reduction does not alter this contribution much. Hence, the preferred viscosity profile is Figure 10. χ 2 misfit for different upper (ν um ) and lower (ν lm ) mantle viscosities using ICE-6G (Argus et al., 2014;Peltier et al., 2015). Results are shown for a lithospheric thickness of 80 km (a), 115 km (b), and 150 km (c). Models containing lower mantle viscosities >10 22 Pa s perform better for all ice histories. The models that have a χ 2 lower than 40 are denoted by circles. Figure 9. χ 2 misfit for different upper (ν um ) and lower (ν lm ) mantle viscosities. Each subplot is made using a different ice history: (a) ICE-6G (Argus et al., 2014;Peltier et al., 2015), (b) LW-6 (Lambeck et al., 2017), and (c, d) two GLAC-1D ice histories (Tarasov et al., 2012). Models containing lower mantle viscosities >10 22 Pa s perform better for all ice histories.
The models that fit the data within 2.5X the standard deviation are denoted by circles.

10.1029/2020JB020484
Journal of Geophysical Research: Solid Earth not sensitive to the conversion factor, in agreement with findings by King (1995), justifying our choice for a single value of the conversion factor of 0.15. Figure 12 shows the gravity field from the best-fitting model and the residual between this model and the gravity observations. We find the lowest misfit for the LW-6 ice history, using an upper mantle viscosity of 4 × 10 20 Pa s and a lower mantle viscosity of 2.56 × 10 22 Pa s. Some residuals can be expected over other areas over North America, as the misfit is only calculated over Hudson Bay. Nevertheless, the negative residual to the southwest of Hudson Bay near the Rocky Mountains deserves special attention because it influences the gravity anomaly inside the region bounded by the contour in Figure 1. There are several possible explanations for this anomaly: Figure 6b suggests that the uncertainty in the density profile is the cause, as a clear

10.1029/2020JB020484
Journal of Geophysical Research: Solid Earth uncertainty over the Rocky Mountains due to the density profile is exhibited. Employing an isopycnal crust indeed improves the fit but does not enable the full removal of the anomaly over the Rocky mountains. Other options are changes in the LAB (see Figure 6d) or the effect of lateral viscosity changes in GIA models (e.g., Geruo et al., 2013;Kuchar et al., 2019;Paulson et al., 2005).

Conclusion and Discussion
In this study, we combined dynamic models for GIA and mantle convection with a CLA model and matched the results to the long-wavelength gravity anomaly. The dynamic pressures caused by GIA and surface dynamic topography are implemented in the CLA model to compute the lithospheric density anomalies that are needed for isostasy. We argue that this is a necessary step to be able to derive a consistent forward gravity field model. Uncertainties in the ice history, crustal model, LAB, and conversion factor are found to be small enough in the long-wavelength domain such that a lower bound can be placed on the lower mantle viscosity.
The best-fitting model to the gravity field observations is found when lower mantle viscosities are larger than 10 22 Pa s. Our results do not constrain the upper mantle viscosity, as the better-performing models are present for the full range of upper mantle viscosities (10 20 to 10 21 Pa s) preferred in previous studies (e.g., Paulson et al., 2007;Sasgen et al., 2012;Tamisiea et al., 2007;Wolf et al., 2006).
Previous studies have suggested that the gravity anomaly over Hudson Bay is mainly due to mantle convection (Cathles, 1975;James, 1992;Peltier et al., 1992). Peltier et al. (1992) found that conversion factors (from seismic velocities to densities) in the range of 0.5-1.5 are needed to explain the gravity anomaly by mantle convection, which is large compared to recent estimates, which are in the range 0.1-0.4 (Trampert et al., 2004). Tamisiea et al. (2007) attributed less than 50% (25-45%) of the anomaly to GIA. They estimated the viscosity based on gravity rates but did not check whether the remaining percentage can be explained by mantle convection and did not include crustal and lithospheric density anomalies. Our results show that, at the modeled minimum, at least 60% of the negative anomaly in the gravity field can be attributed to GIA. In previous studies (Métivier et al., 2016;Simons & Hager, 1997), authors also found a preference for a lower mantle viscosity >10 22 Pa s.
Earlier GIA studies (Caron et al., 2017;Nakada & Okuno, 2016) found that two sets of viscosity profiles result in small misfit values, which is classical for GIA models. The first set of well-performing models contains lower mantle viscosities between 10 21 and 10 22 Pa s, whereas the second set has lower mantle viscosities greater than 10 22 Pa s. Solutions containing lower mantle viscosities between 10 21 and 10 22 are derived from data on RSL (Cianetti et al., 2002), GRACE gravity rates (Tamisiea et al., 2007;van der Wal et al., 2008), GPS (van der Wal et al., 2009), or a combination of two of these (Paulson et al., 2007;Zhao, 2013). ICE-6G is based on the VM5a viscosity structure, which has a lower mantle viscosity <10 22 Pa s (Peltier & Drummond, 2008).
In contrast, several studies have found a high viscosity in the lower mantle, for example, by an inversion of GPS, tide level gauges, absolute gravimetry, and sea-level indicators (Wolf et al., 2006) or by inverting for gravity rate observations from GRACE together with present-day ice mass changes in Alaska and Greenland (Sasgen et al., 2012). Steffen et al. (2009) compared GRACE solutions with results of a GIA model adjusted to fit RSL curves and found 2 × 10 22 Pa s for the lower mantle viscosity. Geological evidence for RSL change and the tilting of paleolake shorelines, combined with present-day crustal movement, converged to high lower mantle viscosity models (Lambeck et al., 2017). Métivier et al. (2016) used gravity gradients and concluded that lower mantle viscosities larger than 2 × 10 22 Pa s are preferred. This agrees with the analysis of _ J 2 data, which required viscosities above 5 × 10 22 Pa s in the lower part of the lower mantle (Nakada & Okuno, 2016). Finally, Kuchar et al. (2019) found that an average viscosity of 3 22 Pa s is needed to fit RSL data in 1-D models and that the evolution of the peripheral bulge near the Atlantic and Gulf coast is what requires these high viscosities. While the area investigated here does not include the peripheral bulge near the Atlantic and Gulf coast, results from our model, constrained by the gravity field, exhibit a clear preference for lower mantle viscosities >10 22 Pa s. The lower mantle viscosity affects inferences based on GIA models, such as the distribution of ice volume required to close the sea-level budget at LGM (Lambeck et al., 2014).
In general, mantle convection studies are global studies, employing a more complex radial viscosity profile than used in our study. Nevertheless, the viscosity found in our study is in rough agreement with studies on 10.1029/2020JB020484 Journal of Geophysical Research: Solid Earth slab sinking speeds (Čížková et al., 2012), mantle convection (e.g., Bower et al., 2013;Perry et al., 2003;Steinberger, 2007), or when mantle convection is combined with GIA (Mitrovica & Forte, 2004). In three-layered mantle convection models, one of the important parameters determining the amplitude and shape of the surface dynamic topography is the increase in viscosity between the upper and lower mantle. Since the upper mantle viscosity over North America is found to lie between 10 20 and 10 21 Pa s (e.g., Paulson et al., 2007;Sasgen et al., 2012;Tamisiea et al., 2007;Wolf et al., 2006), lower mantle viscosities >10 22 Pa s require a jump that is likely to be at least a factor of 20 at the boundary between the upper and lower mantle. This is consistent with other mantle convection studies (Rudolph et al., 2015) but deviates somewhat from a study that found a jump of only 10 between the upper and lower mantle viscosity (Liu & Zhong, 2016). All in all, our results are not in conflict with most studies on mantle convection. Our results also prefer an upper mantle viscosity between 10 20 and 10 21 Pa s and consequently support a contrast larger than 20 between the upper and lower mantle viscosity, which favors slower slab sinking speeds ( Van der Meer et al., 2018).
We have developed an approach to combine CLA models with models for GIA and surface dynamic topography using isostasy. In principle, we solve for the average viscosity value of the complete lower mantle, but most of the sensitivity will be toward the upper part of the lower mantle (e.g., Mitrovica & Peltier, 1991a). Our approach can be applied to other regions that experience ongoing large-scale GIA, like Fennoscandia, Alaska, and Antarctica. The SHs were truncated at degree 15, which diminished uncertainties due to the crustal model that were previously found to be large in Scandinavia (Root et al., 2015). If we studied small-scale GIA signals, a larger uncertainty would be introduced by the crustal model. Since mantle convection covers the long wavelengths, this concept could also be useful for regional mantle convection models that aim to constrain viscosity or the conversion factor from seismic velocities to densities.

Appendix A: Sensitivity to Lithospheric Thickness in the Mantle Convection Model
In Figure 9, we have shown that the main results do not exhibit a large sensitivity to the lithospheric thickness in the GIA model. Besides the sensitivity to the lithospheric thickness in the GIA model, it is also insightful to explore the sensitivity to the thickness of the lithosphere in the mantle convection model. There are also only minor differences in the gravity signal coming from the mantle convection model, up to 0.6 mGal for an upper and lower mantle viscosity of 4 × 10 20 and 1.3 × 10 22 Pa s, respectively ( Figure A1).

Appendix C: The Crustal Signal for lmax 20
In our study, we have only looked at the first 15 SHs. For a maximum degree of 20, the range in the gravity signal is of the same magnitude, but there are more regions where the spread is >10 mGal, as compared to Figure Figure 6 but now for a maximum degree of 20. The spread in gravity field is of the same magnitude, but there are more region with uncertainties >10 mGal, especially in the crust.

Data Availability Statement
Codes to run the crust-lithosphere-asthenosphere model and the GIA model are available at https://doi.org/ 10.4121/uuid:54c182ea-3826-4200-83d6-414dad27b289. The SFEC code for modeling mantle convection is available upon request from Zdenek Martinec.