Uncertainties of Glacial Isostatic Adjustment Model Predictions in North America Associated With 3D Structure

We quantify GIA prediction uncertainties of 250 1D and 3D glacial isostatic adjustment (GIA) models through comparisons with deglacial relative sea‐level (RSL) data from North America and rate of vertical land motion ( U̇ ) and gravity rate of change ( Ġ ) from GNSS and GRACE data, respectively. Spatially, the size of the RSL uncertainties varies across North America with the largest from Hudson Bay and near previous ice margins along the northern Atlantic and Pacific coasts, which suggests 3D viscosity structure in the lower mantle and laterally varying lithospheric thickness. Temporally, RSL uncertainties decrease from the Last Glacial Maximum to present except for west of Hudson Bay and the northeastern Pacific coast. The uncertainties of both these regions increase from 30 to 45 m between 15 and 11 ka BP, which may be due to the rapid decrease of surface loading at that time. Present‐day U̇ and Ġ uncertainties are largest in southwestern Hudson Bay with magnitudes of 2.4 mm/year and 0.4 μGal/year, mainly due to the 3D viscosity structure in the lower mantle.


Introduction
Glacial isostatic adjustment (GIA) studies are vital to interpret relative sea-level (RSL) change, present-day rates of vertical land motion ( _ U), and gravity-rate-of-change ( _ G) signals with respect to the growth and decay of Quaternary ice sheets (e.g., Peltier, 2005;Schumacher et al., 2018). Due to the imperfect knowledge of GIA model inputs (i.e., models of the surface loading history and Earth's internal viscoelastic structure), it is desirable to quantify the uncertainties of GIA predictions from an ensemble of GIA models (Caron et al., 2018;Melini & Spada, 2019). GIA uncertainties have previously been derived from (1) 20% of the magnitudes of the signal of interest (e.g., Karegar et al., 2017;Slangen et al., 2012) for time-dependent gravity measurements of the Gravity Recovery And Climate Experiment (GRACE) twin-satellite mission, (2) using (semi)empirical estimation (e.g., Hill et al., 2010;Riva et al., 2009;Simon et al., 2017), and (3) residual misfit analysis (e.g., Peltier et al., 2015).
Recent advances in computational power have enabled GIA uncertainties to be estimated using a large ensemble of randomly generated ice-Earth models (e.g., Caron et al., 2018;Melini & Spada, 2019). However, the Earth component of GIA models has used 1D (laterally homogeneous) viscosity models that neglect the influence of 3D structure. Surface geology (e.g., Kennett & TkalČić, 2008) and seismic tomography (e.g., Bunge & Grand, 2000) show that Earth's material properties are laterally heterogeneous and many studies have revealed that the 3D structure has a significant influence on GIA predictions (e.g., Austermann et al., 2013;Yousefi et al., 2018). Indeed, 3D structure has been invoked as a mechanism whereby GIA models may better fit late Quaternary RSL records (e.g., Clark et al., 2019;Kuchar et al., 2019;Love et al., 2016), while neglecting 3D structure could introduce bias in 1D viscosity inversions (Lau et al., 2018 The GIA prediction uncertainties of RSL, _ U, and _ G can be separated into two types (Melini & Spada, 2019): (1) uncertainties associated with input parameters of the Earth rheology model or surface mass loading history of continental ice sheets (T1) and (2) uncertainties associated with structural differences among GIA models that use different methods to solve the sea-level equation, such as the implementation or neglect of coastline migration (T2). Here, we subdivide the T1 uncertainties into those associated with Earth model only (T1A), with surface loading history only (T1B) and with both Earth model and surface loading history (T1C) (see Text S1 in the supporting information for detailed discussion on GIA uncertainties and the strategy of this paper). In this paper, we use the surface loading history of ICE-6G_C (Argus et al., 2014;Peltier et al., 2015) and restrict our analysis to North America, where quality-controlled RSL data are densely distributed along the coasts and previous analysis has revealed misfits between 1D GIA models and the data (e.g., Roy & Peltier, 2015). Our goal is to isolate the T1A GIA prediction uncertainties for an ensemble of Earth models. The ensemble has 250 Earth models in total, including 150 3D Earth models. The 3D structure contains lateral variations in mantle viscosity, lithospheric thickness, and sublithospheric and asthenospheric properties ( Figure S1).
To calculate the T1A uncertainties, we select a subset of GIA models from the large ensemble based on the ability of model predictions to fit simultaneously: (1) the latest quality-controlled deglacial RSL databases from the Atlantic (Engelhart & Horton, 2012;Vacchi et al., 2018) and Pacific coasts of North America (Engelhart et al., 2015) and (2) observed present-day _ U and _ G data in North America from Global Navigation Satellite System (GNSS) and GRACE data. Present-day horizontal motion data are not considered here because they have not been used to constrain the ICE-6G_C model . We investigate the GIA prediction uncertainties using 3D global structure, thereby laying the groundwork for future studies that also consider uncertainty in surface loading history (e.g., T1C) and the separation of uncertainties due to surface loading history and Earth model parameters.

Methods
The GIA response of a spherical, self-gravitating, materially compressible Maxwell Earth is computed using the Coupled Laplace-Finite Element method (Wu, 2004). The effects of rotational feedback and time-dependent coastlines are also included in the computation of the solution to the sea-level equation (Clark et al., 1978). Our finite element grid has a 0.5× 0.5°spatial resolution near the surface but decreases with depth to 2.0× 2.0°in the lower mantle to reduce computation time. Note that when the Coupled Laplace-Finite Element method is employed to reconstruct the results for the ICE-6G_C (VM5a) model of Peltier et al. (2015) that was computed with conventional spectral-normal mode method, Text S2 and Figure S2 show that RSL and present-day _ U predictions are very closely recovered.
We search the parameter space of models with varying (1) laterally heterogeneous mantle viscosity, (2) lateral thickness of an elastic lithosphere, and (3) viscosity and lateral thickness of the sublithosphere and asthenosphere that can fit deglacial RSL data and the observed present-day _ U and _ G measurements in North America simultaneously ) (see Texts S3 and S4 and Table S1 for further details concerning the search algorithm). The lateral viscosity variations in the mantle are derived from shear wave velocity anomalies in a global seismic tomography model (Karato, 2008;Wu et al., 2013). We employ the Bunge and Grand's (2000) seismic tomography model with scaling factors to search for 3D viscosity model solutions. We note that use of a different seismic tomography model, S20A of Ekström and Dziewonski (1998), does not significantly affect which combination of parameters and scaling factors produces the best-fitting model . The lateral variation of lithospheric thickness is, however, derived from model S20A (Ekström & Dziewonski, 1998). The importance of lateral thickness variations in the lithosphere (or lack thereof) has been discussed by Nield et al. (2018), Spada et al. (2006), Wu et al. (2010), and others.
We compare GIA predictions to a quality-controlled deglacial RSL database from the Atlantic (Engelhart & Horton, 2012;Vacchi et al., 2018) and Pacific coasts of North America (Engelhart et al., 2015), which is also updated with relevant new published reconstructions (e.g., Hawkes et al., 2016;Khan et al., 2017;Love et al., 2016; Figure S3, Text S4). This combined data set contains 3,368 data points in total, including 1,725 sea-level index points, which define the elevation of RSL at a given point in space and time, and 874 marine limiting and 769 terrestrial limiting data points, which provide lower and upper limits on the position of RSL ( Figure S3), respectively (Engelhart & Horton, 2012;Shennan & Horton, 2002). We selected 24 regions from the Atlantic and Pacific coasts of North America, prioritizing RSL records that are temporally complete across during the Holocene (12,000 years to present, Figure 2a).
We compare present-day _ U and _ G predictions from the ensemble of GIA models with observations processed from GNSS ( Figure 3a) and GRACE data that were used to constrain (GPS observations) and/or validate (GRACE observations) the ICE-6G_C model in Peltier et al. (2015). Additionally, we use the GIA induced _ G signal transformed from GNSS data (see Wang et al. (2013) for methodology, Figure S4) for comparison with the present-day _ G predictions ( Figure 3c).
We calculate the misfit χ statistics between deglacial RSL data, observed present-day _ U and _ G data and the values predicted from 250 GIA models (including 150 3D GIA models) to select the subset of GIA models for uncertainty calculation (Text S4). Among the subset of GIA models, the mean and standard deviation of model predictions of RSL, _ U, and _ G are calculated as follows: ( 2) where N represents the number of best-fitting model selected, p i (m j ) indicates the ith prediction (e.g., at time t i at a given RSL site) of model m j , p i represents the mean of the ith prediction, and σ i represents the standard deviation of the ith prediction. Here, we take 2σ i as the ith prediction uncertainty. Note that our statistics (equations 1 and 2) are different from those previously employed (e.g., Melini & Spada, 2019) and, therefore, the results are not directly comparable.

Results and Discussion
Based on the goodness of fit with the deglacial RSL data and the observed present-day _ U and _ G data in North America, we select the 10 best-fitting models (Text S4 and Table S2) from 250 GIA models (including 150 3D GIA models) to compute the mean and 2σ uncertainties of GIA predictions of RSL, _ U, and _ G. A subset of 5 or 20 best-fitting models does not significantly alter the mean and 2σ uncertainties (Text S4, Figures 1, 3, and S5-S8). The uncertainties of 10 best-fitting models encompass a reasonable range of GIA predictions of RSL, _ U, and _ G with different 3D parameters in the Earth model. The 10 best-fitting models include five 3D models, two 1D models, and three models with 1D mantle viscosity but laterally varying lithosphere, sublithosphere and asthenosphere (Text S4 and Table S2).

Deglacial RSL Uncertainties
The mean of the 10 best-fitting GIA models shows two maxima of RSL in western and eastern Hudson Bay which decrease from the Last Glacial Maximum to present: >600 m at 15 ka BP, >300 m at 9 ka BP, and > 120 m at 6 ka BP (Figure 1). During deglaciation, mean RSL gradually migrates northeastwards from the west of Hudson Bay between 15 and 9 ka BP. This pattern is mainly due to the melting of two ice domes in western and eastern Hudson Bay in ICE-6G_C .
The 2σ uncertainties of the RSL predictions show a different pattern to the mean (Figures 1 and S5) with larger uncertainties west and east of Hudson Bay (~30 m at 15 ka BP) and near previous ice margins of the Laurentide ice sheet along the northern Atlantic and Pacific coasts (~20 m at 15 ka BP). Temporally, the magnitude of RSL 2σ uncertainties decreases from Last Glacial Maximum to present in most locations. However, there is an increase in the magnitude of the RSL 2σ uncertainties from 30 to 45 m between 15 and 11 ka BP ( Figure S5) west of Hudson Bay and the north Pacific coast, which is related to the fast release of surface loading ( Figure S9) from the Meltwater Pulses 1A and 1B (Fairbanks, 1989) in the ICE-6G_C reconstruction that are primarily sourced from the Laurentide ice sheet . Spatially, the patterns of RSL 2σ uncertainties are similar from 15 to 8 ka BP, after which areas of higher uncertainties (15 m at 7 ka BP and ≥1 m at 1 ka BP in Figure S5) emerge from southwest Hudson Bay to southern Saskatchewan and along northern Atlantic and Pacific coasts of North America. Because all ice has melted from North America by~8 ka BP , we suggest these uncertainties may be entirely due to the uncertainties in the 3D viscosity structures in the mantle. Indeed, Love et al. (2016) and Yousefi et al. (2018) concluded that northern Atlantic and Pacific coasts are sensitive to 3D structure. We compared the mean RSL predictions of the 10 best-fitting models and the RSL predictions of the best-fitting model HetM-d110t20v22 (Text S4 and Table S2) with the deglacial RSL data from the 24 selected regions (Figures 2 and S10). The best-fitting model HetM-d110t20v22 includes 3D viscosity structure in the mantle but a lithosphere of uniform thickness-this is because we searched the parameter space of the 3D viscosity structures with a uniformly thick lithosphere . Of the 250 models considered, the model that least fit the data is a 3D model with VM5a as 1D background viscosity model and scaling factors equal 1 ( Figure S11).
The mean RSL with 3σ uncertainties and the predictions of model HetM-d110t20v22 fit most of the deglacial RSL data from near (e.g., regions 1-6 and 11-13) and intermediate (e.g., regions 15-17, 20-21, and 24) study regions (Figures 2b and S10), which validates the performance of our selected models. The RSL uncertainties are relatively small in the near field around Hudson Bay and James Bay (e.g., regions 3-6) compared to the large RSL magnitudes in these regions (e.g., >180 m at 8 ka BP for regions 3-6). RSL uncertainties are also small in the intermediate field distant from the Laurentide ice sheet margins (e.g., regions 17-18 and 22-24, Figures 2b and S10).
The largest RSL uncertainties are found in regions close to the ice margins (e.g. regions 7-9 and 13-15, Figure S10). The locations of these large uncertainties may be related to forebulge migration toward the ice center as the regions experienced both GIA-related uplift early in deglaciation but forebulge collapse later (e.g., Barnhardt et al., 1995;Roy & Peltier, 2015;Stea et al., 2001;Tushingham & Peltier, 1991). Wang and Wu (2006) suggested that forebulge areas (e.g., regions 13-15) are sensitive to the 3D Earth models.  Peltier et al. (2015). The former indicates that our mean present-day _ U with 2σ uncertainties fits with the observed GNSS data at that site, while the latter fits within two times the errorbars. The purple lines (c and d) represent the 54°N cross section shown in Figure 4.

Present-day _ U and _ G Signals
The mean of the 10 best-fitting GIA models show two maxima of present-day _ U and _ G predictions in western and eastern Hudson Bay with magnitudes over 12 mm/year and 1.8 μGal/year (Figures 3a and 3c). The mean _ U and _ G have a similar distribution pattern, decreasing outwards from the two maxima and close to 0 near the ice margins.
The 2σ uncertainties of _ U and _ G are greatest in southwestern Hudson Bay (>2.4 mm/year and 0.4 μGal/year, Figure 3), which are mainly due to lateral heterogeneity in the lower mantle and consistent with recent sensitivity results of . However, the distribution patterns are also affected by the lateral variations of lithospheric thickness, asthenospheric thickness, and viscosity ( Figure S12). For example, the 2σ uncertainties centered offshore along the Pacific coast of central North America are mainly due to a low-viscosity asthenosphere in that area ( Figure S12). The magnitudes of 2σ uncertainties of _ U and _ G are much smaller, and their patterns are different from the results of Caron et al. (2018) (Figure S13) because they applied a wide range of 1D mantle viscosities and also modified the ice history. As a consequence, not all their model predictions can simultaneously fit deglacial RSL histories and observed present-day _ U and _ G data. If we apply the 20% rule of GRACE ( Figure S14), the _ G 2σ uncertainty has a similar magnitude but a different distribution pattern. The 20% rule overestimates the uncertainty in western and eastern Hudson Bay but underestimates the uncertainty near the ice margins, especially in southwestern Hudson Bay where we have the largest _ G uncertainty. Most of the _ G predictions from the 10 best-fitting models along 54°N latitude fall within 2σ uncertainties from the mean _ G (Figure 4). The largest uncertainties extend from 110°W to 85°W which corresponds to the 2σ uncertainty maximum in southwestern Hudson Bay (Figure 3d). Moving east, the trough in _ G near −80°W is located at James Bay and is followed by increasing _ G that peaks around −75°W. The peak _ G predictions are between~1.8 and 2.0 μGal/year and correspond with the peak GRACE observations . The _ G uncertainties then decrease sharply from 70°W to the Atlantic Ocean.
We compared the mean present-day _ U and _ G predictions including their 2σ uncertainties from the 10 best-fitting models (Figure 3) with observations Sella et al., 2007). 339 out of 348 (97.4%) GNSS site observations from Peltier et al. (2015) are fit by the mean _ U within one observed error bar (Figure 3a). The mean present-day _ G (Figure 3c) prediction accurately fit the double "bulls-eye" structure that evident in GRACE time-dependent gravity observations. Also, the mean present-day _ G prediction fits, both in magnitude and pattern, the GIA-induced _ G signal transformed from GNSS data using the method of Wang et al. (2013) (Figure S4).

Conclusions
We provided mean and uncertainties of GIA predictions in North America associated with 3D mantle viscosity profiles (T1A) from an ensemble of 250 Earth models, including 150 3D Earth models that consider lateral variations in mantle viscosity, lithospheric thickness, and sublithospheric and asthenospheric properties. Based on the goodness of fit with the combined RSL data set and the observed present-day _ U and _ G data in North America, we select 10 best-fitting GIA models to compute the mean and 2σ uncertainties of GIA predictions of RSL, _ U, and _ G. The GIA models show the following: 1. The magnitude of RSL 2σ uncertainties varies across North America with larger uncertainties around Hudson Bay and near previous ice margins of the Laurentide ice sheet along the northern Atlantic and Pacific coasts. The magnitude of uncertainties decreases from the Last Glacial Maximum to present except for the increase in magnitude of the maxima off the west of Hudson Bay and northeastern Pacific coast from 30 to 45 m between 15 and 11 ka BP. The uncertainties patterns do not change temporally until 8 ka BP, after which an area of higher uncertainties emerges in southwestern Hudson Bay that may be entirely due to the uncertainty of the 3D structures. 2. The 2σ uncertainties of the present-day _ U and _ G are centered around the southwestern Hudson Bay with magnitudes around 2.4 mm/year and 0.4 μGal/year, which are different from previous uncertainty study results obtained with 1D Earth model-based analyses. 3. The uncertainties of RSL and present-day _ U and _ G strongly imply that the 3D structures need to be considered in studies of deglacial RSL reconstructions (e.g., Engelhart et al., 2009), future sea-level rise projections (e.g., Love et al., 2016), analysis of vertical land motion from tectonics (e.g., Kreemer et al., 2014;Majewski et al., 2018), and signals of hydrology and ice-mass balance changes from GRACE data (e.g., Shepherd et al., 2018;Wang et al., 2013).
In this study, we have only investigated GIA model prediction uncertainties in North America that is associated with the presence of 3D structure (T1A) but with fixed ICE-6G_C surface loading history. Hence, the uncertainties in surface loading history and correlation between loading history and mantle viscosity are not considered, but there is a need for further investigation (Text S1). In the absence of this analysis, the necessity of including 3D internal viscoelastic structure to provide acceptable fits to observational data remains undetermined. Moreover, an updated RSL data set with more extensive spatial and temporal coverage (e.g., Khan et al., 2019) will be required to further improve constraints on the GIA model. All these caveats suggest the necessity of continuing to update the GIA uncertainty analysis associated with presence of 3D structure as new data become available.