Uncertainty Estimation in Regional Models of Long‐Term GIA Uplift and Sea Level Change: An Overview

This work provides a comparison of four approaches that can be used to describe uncertainty in models of the long‐term glacial isostatic adjustment (GIA) process. The four methods range from pessimistic to optimistic representations of GIA uncertainty. Each estimation method is applied to selected one dimensional GIA model predictions and compared with vertical land motion data from Global Positioning System (GPS) measurements across Fennoscandia and North America. The methods are evaluated relative to two main properties: (1) their expected ability to separate non‐GIA from GIA signals and (2) their estimated statistical appropriateness given a specific GIA model and data set. For the first point, non‐GIA signals are considered isolated from the long‐term (millennial time scale) GIA signal at sites where measurement and model uncertainties do not overlap. Across methods, the frequency and accuracy with which non‐GIA signals are separated from GIA signals in GPS data display both consistent similarities and disparities. For the second point, we compare model predictions with rates of vertical land motion and relative sea level change that have been cleaned of non‐GIA signals to determine the most appropriate value of model uncertainty and relate the findings to the four approaches. Best fit inferences suggest that within deglaciation centers, GIA model uncertainty is up to ~2 mm/yr (vertical land motion). Likewise, away from the former ice sheet centers, GIA uncertainty for relative sea level change is inferred to be ~0.3–0.5 mm/yr along the U.S. East Coast and ~0.6–0.8 mm/yr in the North Sea.


Introduction
A well-known issue in modeling of the longer-term viscous glacial isostatic adjustment (GIA) response of the solid Earth is the coupled dependence of the solution on both Earth model parameters and ice sheet parameters. The coupled dependence of GIA predictions on the selected Earth-ice model combination means that model predictions are nonunique; that is, an equally good fit to observations may be achieved for different model combinations (e.g., Schmidt et al., 2014). The nonuniqueness problem can be difficult to eliminate because Earth's rheological structure as well as paleo-ice sheet volume and extent are both indirectly inferred and are therefore only as good as available constraints (the quantity and quality of which vary considerably).
Given the challenge of constraining model input parameters, the uncertainty associated with the GIA process is mostly of a systematic nature, hence difficult to quantify and usually not provided to the users. It is therefore the goal of this study to assess the potential impacts of GIA model uncertainty on the GIA signal at present day. We examine the effect of within-model uncertainty (parameter variation within a model), across-model uncertainty (different GIA model formulations), and methodological uncertainty (how the uncertainty is itself calculated). We distinguish here between long-term or paleo-GIA and contemporary isostatic adjustments driven by recent cryospheric change. The former represents the millennial time scale relaxation response of the Earth to loading and unloading during Pleistocene glaciations (e.g., Peltier & Andrews, 1976;Tamisiea, 2011;Wu & Peltier, 1982), whereas the latter describes the Earth's annual to decadal response to present-day melting of glaciers and ice sheets (Mitrovica et al., 2001). Both long-term GIA and shorter-term isostatic adjustments can be important contributors to present-day vertical motion and sea level change at local scales.
Increasingly, researchers seek to understand the impacts of a warming climate and anthropogenic forcing on the cryosphere and oceans, the hydrological cycle, and coastal sea levels (e.g., Jevrejeva et al., 2018;Larour et al., 2017;Slangen et al., 2012;Wada et al., 2010). Climate impact assessments are impeded by incomplete ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. and uncertain knowledge of climate processes, and several narrower goals emerge, among them, the need to separate accurately signals arising from paleo-GIA from those caused by various contemporary processes. The Gravity Recovery and Climate Experiment (GRACE, Tapley et al., 2004) in particular has driven this need: GRACE has an unprecedented ability to track Earth's surface mass redistributions, but the accuracy at which it can do so is reliant on predictions of paleo-GIA that correct for the effect of Earth's internal mass redistribution (Ivins et al., 2013;Shepherd et al., 2012;Velicogna & Wahr, 2006;Wang et al., 2013;Wouters et al., 2014). This limitation is not restricted to GRACE data; it rather holds for numerous data sets (e.g., Davis & Mitrovica, 1996;Riva et al., 2017) and can affect interpretations of cryospheric change, sea level variation, and terrestrial water storage. However, despite a general interest, studies rarely discuss in detail the uncertainty of GIA model predictions (e.g., Caron et al., 2018;Melini & Spada, 2019) and seldom examine the implications of using different estimations of model uncertainty. It is therefore important to consider how best to define and quantify GIA model uncertainties as the estimates may impact a variety of applications, including present-day sea level change and ice sheet melting.
In this paper, two broad aspects of GIA uncertainty are discussed. First, the GIA uncertainty associated with the nonuniqueness of GIA models (section 2.1) is considered through the comparison of different published models. Second, the method by which GIA uncertainty itself is represented is examined-for this, four methods, varying from pessimistic to optimistic and from formally statistical to back of the envelope, that researchers have used to consider long-term GIA uncertainty (section 2.3) are compared. Three main sections are dedicated to illustrating the impact of GIA uncertainties in scenarios relevant to GIA research applications: The ability to separate GIA from non-GIA signals (a common application of GIA model users) is assessed (section 3), and the viability of each uncertainty estimation method is evaluated within the context of vertical land motion (VLM) rates (section 4) and relative sea level (RSL) change rates (section 5).
The time window of interest is present day; thus, the focus here is the ability to separate the paleo-GIA signal occurring at present day from remaining contemporary signals (or conversely, to identify regions where residual signals may be masked by model uncertainty). We use VLM data, which provide high-precision pointwise measurements of surface deformation as well as RSL change rates. The study areas are Scandinavia and North America, although uncertainty in long-term GIA estimates will also influence mass loss estimates in Antarctica and possibly Greenland (Khan et al., 2016;King et al., 2012;the IMBIE team, 2018;van Dam et al., 2017;Wake et al., 2016). Focusing the tests outside of Greenland or Antarctica allows the sensitivity of the various uncertainty estimation methods to be assessed for a broader range of deformation behavior, as VLM varies from strongly GIA-dominated to strongly dominated by other processes. The goal of the paper is to examine the impact of GIA model uncertainties within the context of common research applications, not to suggest a preferred GIA model out of those selected.

Description of Models, Data, and Uncertainty Methods
This section describes the main components of the study. Present-day predictions of three GIA models are used (section 2.1). As well, four possible methods of representing model uncertainty are discussed (section 2.3). The impacts of varying the model, as well as the estimation method, are examined by applying the uncertainties to the predictions-this yields several combinations of model predictions and model uncertainty, which can be compared with measured signals. Vertical land motion rates are the observational constraints of the first two applications (sections 3 and 4) and are described in section 2.2. The selected models represent a subset of available GIA models, and the uncertainty methods can differ outside of the variations considered here (either due to regional data availability or internal parameter variation); thus, the results are limited by the choice of model and method. The results of course also depend on the data set used for comparison: Different observables, longer time series, and regional coverage (which can be limited in Antarctica and Greenland in particular) may all influence the results. For the time being, we view the selected models, data, and methods to provide an adequate cross section of scenarios.

Models
No model can perfectly capture the true GIA signal, but all GIA modeling studies are based on the premise that glacial isostatic deformation can be estimated in a useful way by a model. Thus, the thought experiment begins with the assumption that the vertical GIA signal is known to within some uncertainty bounds; the goal is then to determine the consequences of applying the uncertainty bounds of different uncertainty estimation techniques to the modeled signal. We examine three variants of the modeled GIA signal over Scandinavia and two variants over North America: ICE-6G_C (Argus et al., 2014;Peltier et al., 2015), a variation of the Australian National University (ANU) ice sheet model (considered only for Fennoscandia) (Lambeck, 1995;Lambeck et al., 2010), and the D1 semiempirical model (Simon et al., , 2018 (Figures 1a-1e). Vertical uplift predictions for ICE-6G_C are taken directly from online (http://www. atmosp.physics.utoronto.ca/~peltier/data.php) and are coupled to the VM5a rheological profile which has a 60 km thick elastic lithosphere underlain by a 40 km thick high viscosity (10 22 Pa s) layer, which are both then underlain by 5 × 10 20 Pa s upper mantle and a lower mantle that transitions from approximately 1.6 × 10 21 to 3.2 × 10 21 Pa s viscosity below around 1,200 km depth. Relative to ICE-6G_C, VLM predictions from ICE-6G_D can vary significantly in parts of Antarctica (Peltier et al., 2018;Purcell et al., 2016); however, the predictions of the two model versions do not vary significantly over the study areas, so the results will not be impacted by the selected model version. The second ice sheet reconstruction uses the glacial history of the ANU model(s) over Scandinavia and the British Isles; elsewhere on the globe, the model is filled in with ICE-5G (Peltier, 2004) and the North American ice sheet coverage from Simon et al. (2015Simon et al. ( , 2016. Tests indicate that predicted regional rates are not particularly sensitive to the version(s) of the ice sheet reconstruction used outside the study area (the average difference in predicted uplift rate is <0.05 mm/yr magnitude over the study area). Predictions for the ANU reconstruction were obtained by pairing the ice sheet model with Earth model parameters indicated by Lambeck et al. (2010) to be within the range appropriate for Fennoscandia, specifically, a 90 km thick elastic lithosphere and respective upper and lower mantle viscosities of 4 × 10 20 and 1 × 10 22 Pa s. The D1 model is the result of inverting GPS-measured crustal velocities simultaneously with a suite of Earth ice model combinations (Simon et al., , 2018 and thus has no fixed description of Earth structure underlying it. Note than none of the selected models considers changing ice coverage in the Alps. The predicted patterns of uplift are broadly similar for the three models (Figures 1a-1e). In Scandinavia, all models predict peak uplift rates in the Gulf of Bothnia of up to 8 mm/yr, although the spatial location of peak uplift differs slightly. The models also indicate uplift in the Barents Sea, although ICE-6G predicts higher uplift rates there than either of the ANU or D1 models. ICE-6G and D1 both predict uplift over the northern British Isles, whereas in the ANU model, the uplift associated with the British Isles Ice Sheet is masked by the subsidence caused by the peripheral forebulge of the Fennoscandian Ice Sheet. It is worth noting that the ANU model assumes a thinner lithospheric thickness in the British Isles (e.g., Lambeck, 1995); incorporation of varying lithospheric thickness may yield somewhat different predictions than shown in Figure 1. In North America, ICE-6G predicts a three-dome pattern of uplift east and west of Hudson Bay and over Foxe Basin. Compared to ICE-6G, D1 predicts similar peak uplift rates east of Hudson Bay, but without a prominent three-dome pattern.

Vertical Velocities
In this study, vertical velocities are obtained from Global Navigation Satellite System (GNSS) network measurements, specifically from the Global Positioning System (GPS). In Scandinavia, rates of VLM measured by GPS are a combination of data from both Kierulf et al. (2014) and the Nevada Geodetic Laboratory (NGL) (Blewitt et al., 2016). The denser coverage within the former load center in the Kierulf et al. (2014) data set is elsewhere supplemented with data from Blewitt et al. (2016) (Figure 1f).
In North America, rates of VLM are taken from the Blewitt et al. (2016) NGL data set, and are filtered for length, data quality, and spatial overlap in a similar manner as the Fennoscandian rates ( Figure 1g). Extended description of the GPS data is provided in the supporting information.
The GPS measurements are also grouped into three conceptual categories: GIA dominant, non-GIA dominant, and mixed ( Figure 2 and Figures S1-S3 in the supporting information). The non-GIA VLM signal is inferred from the measured signal and modeled GIA rate and the magnitude of the ratio of non-GIA:GIA VLM determines the classification. Those sites at which the ratio of non-GIA:GIA VLM is <0.5 are classified as GIA dominant. These sites are usually located within the former load centers of the Laurentide and Fennoscandian ice sheets and away from Greenland and present-day glaciers; it is the expectation here that VLM rates recorded by GPS will be dominated by the paleo-GIA signal. Sites at which the ratio of non-GIA: GIA VLM is >1.5 are classified as non-GIA dominant. It is the expectation that VLM at these sites is dominated by processes other than long-term GIA, such as tectonics, contemporary GIA, groundwater or hydrocarbon withdrawal, natural hydrological variations, and/or sediment loading. The remainder of the sites have non-GIA:GIA VLM ratios of 0.5-1.5 and are classified as mixed. These sites may contain an unknown mixture of long-term GIA, contemporary GIA, and non-GIA signals. The bounds of the classification scheme are essentially arbitrary, but the distinction provides a convenient way to view broad trends across the data set. In sections 3 and 4, all GPS measurements are compared in the same way, regardless of classification (i.e., the distinction is visual, not computational). The classification is applied separately to the three GIA models, leading to slightly different classifications for each, although prominent GIA dominance in the former Scandinavian and North American load centers remains a constant feature ( Figure 2).

Model Uncertainty Estimation Methods
We examine four methods of uncertainty estimation that have been applied or discussed in various GIA studies (Table 1): (i) parameter variation, (ii) residual analysis, (iii) the canonical ±20% value, and (iv) (semi) empirical estimation. There are likely additional approaches, but these four methods span what can be considered as likely pessimistic and optimistic representations of GIA uncertainty. Parameter variation considers various Earth ice model combinations, the standard deviation of which can be interpreted as a measure of uncertainty (Caron et al., 2018;Wu et al., 2010); this method captures a wide range of possible GIA behaviors but may be overly pessimistic, particularly in load centers. Moreover, there are multiple combinations of parameters that could be considered for variation; selection of which parameter(s) to vary and by how much is itself a source of uncertainty. Residual analysis considers the fit of the GIA model predictions to a set of constraining data. This method does not constitute a true measure of uncertainty, but it is commonly used to assess the ability of a given model to explain long-term GIA (Peltier et al., 2015). The third method is a rule of thumb, which assumes that in general the uncertainty associated with GIA is within approximately ±20% of the total GIA signal. This method possibly originated in the analysis of Paulson et al. (2007) who presumed a rather heuristic ±20% uncertainty in models of ice sheet thickness around Hudson Bay; the authors allowed for ±20% ice sheet model uncertainty by scaling predicted model outputs to vary within Purple circles-GIA dominated, green diamonds-non-GIA dominated, and orange triangles-mixed. As described in the text, the classification scheme is defined separately for each GIA model. these bounds. The ±20% quantity has since sometimes been adopted when applying uncertainty to the GIA correction for GRACE measurements or rates of sea level change (e.g., https://grace.jpl.nasa.gov/data/getdata/gia-trends/, Karegar et al., 2017;Slangen et al., 2012). The final method is semiempirical estimation, which generally involves the inversion of constraining data sets to arrive at a formal measure of GIA uncertainty (Hill et al., 2010;Riva et al., 2009;Sasgen et al., 2012;Simon et al., 2017).
Since forward models generally do not have a formal measure of uncertainty provided, only the first three methods are applied to the predictions from the ICE-6G and ANU models (Table 1). Likewise, because the D1 model predictions already have a suite of GIA models and the GPS measurements incorporated a priori, neither the parameter variation nor the residual analysis methods are applied to the semiempirical D1 estimates (Table 1). ICE-6G is likewise based on GPS measurements, but we still perform the residual analysis for this model because the model also contains independent constraints and the GPS time series are now longer than when ICE-6G was developed.
Figures 3 and 4 show uncertainty maps for the four estimation methods. Uncertainty values for the parameter variation are obtained by applying the relevant ice sheet history (either ICE-6G or the ANU model) to a set of Earth models that fixes the elastic lithospheric thickness at 90 km and varies upper mantle viscosity values from 0.2-2 × 10 21 Pa s and lower mantle viscosities from 1-60 × 10 21 Pa s; the standard deviation of the model set provides the uncertainty estimation. While varying the viscosity profile relative to a fixed ice sheet reconstruction may violate the near and far-field sea level constraints used to tune the ice sheet reconstruction, the ice sheet model itself can be considered as a single parameter in a larger set of possible models; for practical purposes, what we seek here is a snapshot of possible uncertainties generated by parameter variation that indicates large (pessimistic) uncertainties that scale with the GIA signal (Figures 3 and 4). The parameter variation differs somewhat between the ICE-6G and ANU models due to differences in their loading history, although the spatial pattern of uncertainties is broadly similar. For the residual analysis, the residual between the model predictions (Figures 1a-1e) and data (Figures 1f and 1g) is calculated at each site (Figures 3 and 4) and then the weighted average residual is applied as the uncertainty; unlike the other estimation methods, this representation of uncertainty is spatially invariant. Non-GIA-dominated and mixed sites are excluded from the calculation of the average residual, since there is no expectation that the model should fit the data well at sites not dominated by GIA. The ±20% uncertainty values are obtained as the magnitude of the model predictions in Figures 1a-1f multiplied by a factor of 0.2. The semiempirical uncertainties are from Simon et al. (2017Simon et al. ( , 2018 and compare well in magnitude to similar semiempirical estimates from Hill et al. (2010), as well as those presented as part of the Stable North American Reference Frame (SNARF) project (https://www.unavco.org/projects/past-projects/snarf/snarf.html). In the former load May be unnecessarily large, particularly in load centers (e.g., Table 2); selection of which parameter(s) to vary and by how much itself subject to uncertainty 2. Residual analysis Provides wholesale assessment of quality of 1 GIA model relative to body of data Not true quantification of GIA uncertainty since only 1 modeled signal is considered 3. Within ±20% of GIA signal Straightforward, few assumptions, accessible also to non-GIA modelers Possibly too large in load centers and too small in far field, will approach zero at hinge line 4. (Semi)empirical estimation Combines advantages of 1) and 2); relatively or completely free of forward modeled GIA uncertainty Still dependent on input assumptions and available data, perhaps optimistic in load centers

Summary of comparisons performed in this study
Journal of Geophysical Research: Solid Earth centers, the semiempirical uncertainties are up to an order of magnitude smaller than those of the other estimation methods. In Scandinavia, the spatial pattern of the semiempirical uncertainties differs from the parameter variation and ±20% methods because the semiempirical uncertainties are largest not only where the signal is largest but also where the data are sparsest and/or more poorly constrained.

Application 1: Separating GIA From Non-GIA Signals
Once obtained, the different model uncertainty estimations are applied to each of the modeled GIA signals following Table 1 and compared to the GPS-measured vertical rates (which have independent uncertainties). GIA and non-GIA signals are deemed resolvable from each other in the data if the measured total rate and the modeled GIA rate do not overlap within their respective 1σ uncertainties. The difference represents the minimum resolvable value of the non-GIA signal (that is, the resolvable non-GIA signal could be larger, and the minimum is a lower bound Figures 5 and 6). We also compute the average resolvable non-GIA signal ( Figures S4 and S5). In this application, there is no assumption of a known non-GIA signal that we aim to quantify correctly; rather, we wish to demonstrate how the method of GIA uncertainty selection impacts the ability to identify non-GIA signals, should they be present. This application holds the data uncertainties fixed and varies the model uncertainties; longer time series and smaller data uncertainties should improve the ability to separate non-GIA and GIA signals.

Fennoscandia
The ability of the various methods and models to resolve non-GIA signals from GIA signals is summarized in Figures 5, 6, and S6 and Table S1, in terms of both the conceptual classification and the minimum value of the resolvable non-GIA signal. For all methods and models, GIA and non-GIA signals can be resolved from each other at between 32% and 56% of the selected GPS sites. As expected, non-GIA signals are resolved least often at GIA-dominated sites (usually <10%) and most often at non-GIA-dominated sites (≥68%) for all model  combinations (Table S1); this tendency is a consequence of regions with larger signals (i.e., former load centers) having larger model uncertainties (and thus less resolving ability) than regions with smaller signals.
In Table S1, the qualitative rankings between models influence the interpretation more than the quantitative percentage values; note that if this were a global study rather than a regional study focused on former Pleistocene ice sheet centers, the total percentage values would change to reflect the inclusion of proportionally more non-GIA sites.
Across methods, the parameter variation resolves the smallest number of non-GIA signals and shows almost no ability to resolve signals at mixed or GIA-dominated sites. The residual analysis resolves a similar amount of non-GIA signals for the ICE-6G and ANU models, although the spatial distribution of the resolved sites differs in the former British Isles and Fennoscandian load centers, reflecting the differences between the two sets of model predictions. Conversely, for all models, the ±20% method shows a consistent spatial distribution of resolvable sites inversely proportional with distance from the load center (although tests in sections 4 and 5 indicate the magnitude of these uncertainties is too small). The semiempirical method performs similarly to the ±20% method for D1, except that the semiempirical method resolves noticeably more non-GIA signals in GIA-dominated regions due to smaller model uncertainties there (Figures 5, S6, and S7). Because the Scandinavian and British Isles component of the ANU ice sheet model have different best fit Earth models, the pattern of resolved non-GIA VLM for our 1-D formulation of the ANU model is probably less reliable over the British Isles than in Scandinavia.
The various scenarios also share common features ( Figure 5). For example, all scenarios resolve non-GIA signals of at least +0.4 mm/yr in southeast England. In the Netherlands, even though predicted GIA rates and measured GPS rates are relatively small, all but two scenarios can resolve negative non-GIA signals at three to four of the four selected sites. All scenarios additionally show an envelope between the zero contour of the GPS data and the zero contour of the GIA model predictions in which non-GIA signals of between −0.4 and −2.0 mm/yr are resolved for between four and six sites. In the southernmost part of the study area, predicted GIA signals are small, and the models resolve a mixture of positive and negative non-GIA signals, the magnitudes of which are reasonably consistent across methods and models (although note that none of the models consider paleo-ice cover in the Alps, which may lead to underpredicted GIA motions in that region). In D1, the semiempirical method resolves positive non-GIA signals at several sites in the former load center with consistent magnitudes of up to +0.4 mm/yr, possibly allowing the detection of the small positive signal expected from elastic uplift from Greenland mass loss .

North America
The results for North America are broadly like those for Fennoscandia (Figures 6, S6, and S8 and Table S2). Perhaps the most noticeable difference is that GIA and non-GIA signals are resolved from each other at between 57% and 76% of the selected GPS sites or at roughly 25% more of the sites than in Fennoscandia (although again the qualitative rankings are of more interest than the exact percentages). The~25% increase likely reflects the fact that there are proportionally more non-GIA sites than GIA-dominated sites in North America compared to Fennoscandia, which again indicates that the total percentages are sensitive to the regional scope of the study area. As with Fennoscandia, non-GIA signals are resolved most often at non-GIAdominated sites (>90%), and least often at GIA-dominated sites (<10%) (Figures S6 and S8 and Table S2).
Again, there are some consistent features among the five scenarios ( Figure 6). All methods resolve significant non-GIA uplift along the northwestern coastline of the continent, with the ±20% and the semiempirical methods resolving slightly larger signals than either the parameter variation or residual analysis. These signals are consistent with uplift from present-day and Little Ice Age mass loss in Alaska and tectonic deformation along the active subduction margin. Significant non-GIA subsidence is resolved in the southern part of the study area around the Gulf of Mexico and at several sites in the southwestern United States, consistent with local sediment loading and/or subsurface fluid withdrawal. These uplift and subsidence trends are so far expected, since the above-mentioned processes likely contribute large non-GIA signals to measured VLM rates in these regions (Harig & Simons, 2016;Ivins et al., 2007;Karegar et al., 2015;Kolker et al., 2011;Larsen et al., 2004;Mazzotti et al., 2007;Montillet et al., 2018;Riel et al., 2018;Wada et al., 2010). Perhaps more interesting are the non-GIA trends that are present at sites along the eastern coastline and within the central continent and former load center. Along the eastern coastline, all scenarios resolve non-GIA signals in VLM at a maximum of four sites. In the central continent south of the Great Lakes, all scenarios resolve a region of non-GIA subsidence; here, subsidence associated with movement of the peripheral bulge in either of the ICE-6G or D1 models appears unable to fully explain the measured GPS motions.

Application 2: Evaluating the Uncertainty Methods With VLM Rates
Section 3 illustrates how much starting assumptions may impact the ability to separate GIA from non-GIA signals but offers little guidance to modelers (or model users) on how best to estimate (or apply) GIA model uncertainties. In this and the following section, the performance of selected uncertainty methods is evaluated quantitatively using a χ 2 analysis.
The normalized χ 2 measure of misfit between N observations and predictions is computed as where for the ith point, v data and v model are the VLM rates of the data and model, respectively, and σ tot represents the total uncertainty. While many forward modeling studies of GIA provide an acknowledgement of GIA uncertainty in the form of parameter variation in Earth ice model combinations (e.g., Caron et al., 2017;Lambeck et al., 1998;Milne et al., 2004;Tamisiea, 2011), most predictions of a single

Journal of Geophysical Research: Solid Earth
GIA model are provided without quantitative uncertainty estimates. However, explicitly including both GIA model uncertainty (σ model ) and data uncertainty (σ data ) yields an expression for total uncertainty σ tot ¼ σ 2 data þσ 2 model À Á 1=2 (similar in approach to Whitehouse et al., 2012, andGunter et al., 2014). An appropriate trade-off between misfit and variance is achieved when χ 2 ≈ 1. A value of χ 2 ≫ 1 indicates that either the uncertainties are underestimated or that the data are not well described by the model, whereas a value of χ 2 ≪ 1 suggests overly conservative uncertainties or an overparameterized model. As in section 3, sections 4 and 5 hold the data uncertainties fixed and vary the model uncertainties and predictions; the χ 2 results therefore represent a view of expected model uncertainty as indicated by recent data and may evolve as data coverage and accuracy improve.
In what follows, we use a subset of both the models and the GPS data introduced in the previous section to compute χ 2 values for Scandinavia and North America. Here, the values of v data and σ data in Equation 1 are taken from Schumacher et al. (2018), who presented a global data set of GPS-measured VLM rates in which the rates have been cleaned of non-GIA signals using a consistent methodology. Schumacher et al. (2018) evaluated the cleaned GIA data with respect to GIA model predictions by calculation of regional root mean square values, which do not include data or model uncertainties. Although the authors are appropriately cautious in noting that these data cannot perfectly represent VLM due to GIA, the Schumacher et al. (2018) data set appears to robustly capture the GIA process. These data therefore provide a convenient input for Equation 1, which, for consistency, requires that both v data and v model describe VLM due to GIA alone. Thus, unlike section 3, which required the full GPS rates to examine the extent to which GIA and non-GIA signals could be separated, the χ 2 analysis utilizes cleaned rates to evaluate directly the viability of the GIA uncertainty methods. The original rates in Schumacher et al. (2018) also come from the NGL database and have been subsequently filtered to remove non-GIA signals. As part of their selection process, the authors exclude locations where VLM likely does not significantly reflect GIA (for example, Gulf of Alaska sites are excluded). Therefore, the χ 2 calculations below use a somewhat different but often overlapping set of the locations shown in sections 2 and 3.
Assuming that the cleaned GIA data are accurate to within their uncertainties, the χ 2 analysis should indicate which representation of GIA uncertainty gives the most statistically appropriate χ 2 value. Using the uncertainty maps described in section 2.3 (Figures 3 and 4) as inputs of σ model in equation 1, we compare the χ 2 values obtained for the parameter variation, the ±20% and the empirical estimation methods. We also examine the case where only the data uncertainties are used in the χ 2 calculation (σ model = 0).
The results for Scandinavia and North America when ICE-6G is the starting GIA model are summarized in Table 2. The corresponding results for the ANU and D1 models are summarized in Tables S3-S5 and Figure S10. Both the regions as a whole and the three regional subdivisions (GIA, mixed, and non-GIA) are considered; the subdivisions are still determined using the uncleaned rates (section 2.2 and Figures 1f  and 1g) since the data from Schumacher et al. (2018) formally speaking should all fall within the GIA Table 2 The χ 2 Results for Scandinavia (Upper)

10.1029/2019JB018983
Journal of Geophysical Research: Solid Earth dominated category. We also limit the calculations for the GIA dominated region to locations with model predicted GIA rates of >0 mm/yr in order to restrict the calculation to the region north of the hinge line. The χ 2 misfit values logically decrease as larger model uncertainties are incorporated (from left to right in Table 2); this trend holds for both the regions as a whole and the three regional subdivisions. Conversely, the χ 2 values consistently increase moving outward from the load centers to the farther-field areas. Three factors contribute to this trend: Moving away from the load center, there are generally (1) lower signal-tonoise ratios, (2) somewhat smaller regionally averaged data uncertainties, and (3) when applied, decreasing model uncertainties. With χ 2 values of ≤0.3 the parameter variation method is the only method where the model uncertainties are ever significantly overestimated and then only within the load centers. Within the load center, the ±20% method gives χ 2 ≤ 2 for both continents, so we search for the percentage that yields χ 2 = 1 ( Figure S9). A χ 2 value of 1 for ICE-6G is achieved at ±32-33% in Scandinavia and ±25-26% in North America (±30% is also shown in Table 2 for comparison). These values are consistent with averaged regional uncertainties within the load center of σ m = ±2.1 mm/yr and σ m = ±1.7 mm/yr for Scandinavia and North America, respectively. The D1 model (Tables S3 and S5 and Figure S10) follows similar trends, except that the empirical estimation method appears too optimistic in both the Scandinavian and North American load centers. For the ±% method, χ 2 values of 1 are achieved for D1 predictions with ±35-36% and ±21% for Scandinavia and North America, respectively. These percentage values correspond to average regional model uncertainties of σ m = 1.7 mm/yr and σ m = 1.0 mm/yr; the latter is smaller than that inferred for ICE-6G, although it is still larger than the corresponding empirical uncertainties.

Application 3: Evaluating the Uncertainty Methods With RSL Rates
In this section, we repeat the test of section 4 using selected present-day rates of RSL change. We focus on the coastlines of the eastern United States and the North Sea. The χ 2 calculation follows Equation 1, where v data and v model are replaced by the RSL equivalents RSL data and RSL model . Frederikse et al. (2016Frederikse et al. ( , 2017 examined the sea level budget at selected tide gauge locations along the U.S. East Coast and in the North Sea (Figure 7). The authors computed present-day mass and ocean dynamics contributions to RSL change, which we remove from the observed total RSL rates to obtain a cleaned, GIA-related RSL signal at these locations (RSL data ). The present-day mass contribution consists of predicted sea level fingerprints for present-day ice melting (from glaciers and ice caps and the Greenland and Antarctic Ice Sheets), as well as modeled changes to continental hydrology. The ocean dynamics contribution along the U.S. coastline and North Sea shelf are approximated by estimates for high-latitude steric variability in the North Atlantic. The estimation of the present-day mass and ocean dynamics contributions is described further in Frederikse et al. (2016Frederikse et al. ( , 2017.

U.S. East Coast
Regardless of which GIA model prediction is applied, the paleo-GIA signal comprises a significant fraction (on average ≥30-50%) of the total observed RSL change at these locations. In Figure 7, the GIA-related RSL rates are compared to model predictions from ICE-6G, D1, and the D3 model (the latter contains a combination of GRACE and GPS data  and is used in the budget calculations of Frederikse et al., 2017, so we include it here as well). Within the empirical uncertainties, the D1 and D3 models are indistinguishable from each other except at Sites 12-14. From north to south along the profile, migration from load adjacent regions to the peripheral forebulge area is reflected by a transition from decreasing RSL rates to increased rates of RSL rise (Figure 7).
Comparing common methods, the ±20% estimation gives χ 2 values that are highest for ICE-6G and lowest for D3. This trend holds for the mean square error (the mean of the sum of the squared residuals) of the three models (Table 3), so it is not simply a result of the ICE-6G model having smaller ±20% uncertainties and thus larger χ 2 values. When the uncertainties from parameter variation are used with the ICE-6G model, a χ 2 value <1 is obtained with an inferred model uncertainty of σ m = 1.0 mm/yr (Table 3). However, these uncertainties are so large that they generally cannot distinguish between any of models ( Figure 7). And while the parameter variation method yields the smallest χ 2 value with ICE-6G, this model also simultaneously has the largest mean square error, suggesting that parameter variation is simply too pessimistic for this region. Conversely, there is little difference between the ±20% method and the empirical estimation 10.1029/2019JB018983 Journal of Geophysical Research: Solid Earth method for the D1 and D3 models, indicating that (semi)empirical type uncertainty estimation may be quite appropriate in this region. Based on either of the D1 or D3 models, it appears that appropriate GIA model uncertainty for this region is σ m~0 .3-0.5 mm/yr (Table 3). Because the average data uncertainty is almost always larger than the inferred model uncertainties, improvement of measurement and/or non-GIA correction uncertainties in this region may also be helpful (Table 3).

North Sea
Rates of GIA-related RSL change are on average smaller in the North Sea than along the eastern U.S. coast (Figure 7) and the modeled GIA rates fit the data somewhat worse, particularly at the more northerly sites. Especially at station Lerwick (Site 10) the estimated contributors may not adequately capture possible influence of the deeper ocean on RSL rates. Regionally, both the ±20% and empirical uncertainty values underestimate GIA model uncertainty. The regional GIA uncertainty is at least equal to the data uncertainty (σ m ≥ σ d ) with calculations with Equation 1 indicating that σ m rates of~0.6-0.8 mm/yr are needed to fit the cleaned RSL data. The larger model uncertainty rates are consistent with the fact that there are still significant across-model variations in GIA predictions in the North Sea (Vermeersen et al., 2018) and indicate that in this region the across-model GIA uncertainty likely dominates the uncertainty of any given model.

Discussion
Using preexisting models and data, we have compared four methods that assess the uncertainty and quality of GIA model predictions. All four methods perform in a consistent manner making them all potentially acceptable candidates for uncertainty estimation. However, the residual analysis is not a true measure of uncertainty so perhaps serves best as a more general evaluation of individual GIA models, particularly since its formulation here is spatially invariant and therefore likely unrealistic.
In section 3, we demonstrate the frequency and accuracy with which non-GIA signals can be expected to be resolved from the paleo-GIA signal, under the assumption of a given method or starting GIA model. For a given model, the number of resolvable non-GIA signals generally increases from parameter variation (Method 1) to semiempirical estimation (Method 4), which is a logical outcome given the reduction in (peak) GIA uncertainty values from Methods 1-4. Across methods, the magnitude of the identified non-GIA signals changes minimally for a given model, although the ±20% and semiempirical methods resolve somewhat larger magnitudes than parameter variation (Figures 5 and 6).
In sections 4 and 5, we use GIA data (i.e., data cleaned of non-GIA signals) to evaluate the performance of the four methods more quantitatively and estimate regionally averaged model uncertainties for regional cases. Generally, the parameter variation method as used here appears overly pessimistic for 1-D models of GIA (Table 2); however, we note there are many possible formulations for parameter variation, which have not been applied in this study. The use of three-dimensional Earth models or different rheological models may also further complicate the expression of GIA uncertainty with the introduction of more free parameters and extrapolation from laboratory to natural scales (e.g., Latychev et al., 2005;van der Wal et al., 2010;Wu et al., 2010). In contrast, the semiempirical method often indicates overly optimistic model uncertainties; future studies could focus on optimizing this approach, as it is the only formal statistical methodology used. In section 5, the RSL analysis of the U.S. East Coast suggests that the semiempirical method currently appears to perform best in farther field regions.
The results of sections 4 and 5 also show the sometimes-quoted nominal GIA uncertainty value of ±20% may be too optimistic in current models of long-term GIA. When compared with Schumacher et al.'s (2018) VLM data within the North American and Scandinavian deglaciation centers, GIA model predictions generally require ±30% (or up to~2 mm/yr) associated model uncertainty to obtain optimal χ 2 values. It seems robust in deglaciation centers that model uncertainties exceed data uncertainties so application of a model uncertainty should not be neglected. The origin of the ±20% method for applications within load centers may mean that this method is not suitable when applied elsewhere, although this method is sometimes used in global reconstructions of sea level.
The regionally averaged GIA uncertainties presented in sections 4 and 5 are for the present-day VLM and RSL signals (~±2 mm/yr VLM for deglaciation centers and < ±1 mm/yr RSL at two locations outside of deglaciation centers). The regional uncertainty estimates are calculated as a sort of empirical "postprediction" estimate and therefore are rather superficial as they do not explicitly describe "paleo" GIA uncertainty and cannot decouple Earth model from ice model uncertainty. The χ 2 evaluation does, however, provide a  . σ d ; σ m ð Þare regionally averaged data and model uncertainties, respectively. Model uncertainty is zero for the "data-only" scenario.

10.1029/2019JB018983
Journal of Geophysical Research: Solid Earth quantitative estimate of the regional present-day GIA uncertainty associated with any given model, as well as point toward a suitable method by which to generate the uncertainty values. In this sense, the estimates are useful since they are directly relatable to what researchers do when they apply a GIA correction for a given model. The inferences are of course dependent on the combination of the data and the model, but they provide a starting basis for smaller regional studies for expressing GIA uncertainties that are appropriate for the selected model and/or region.
Finally, we revisit the starting assumption, which in all cases is that the paleo-GIA signal is known a priori to within some uncertainty. This assumption is implicit whenever a GIA model prediction is utilized; however, it is always possible that the starting premise is incorrect. That is, that resolved non-GIA signals may reflect a deficiency in the starting GIA model, rather than deformation from another process, a fact which may be particularly true in formerly glaciated regions where the GIA signal is larger. Sometimes, the results are sensitive to the selection of starting GIA model, which reflects the limitation of assuming the correctness of any given GIA model. Melini and Spada (2019) recently also indicated significant across-model GIA uncertainties that may effectively be larger than internal GIA uncertainties within a fixed model.
In summary, the impact of using different estimations of GIA uncertainty in one-dimensional GIA models has been thus far seldom examined or quantified and there exists no consensus on how to best represent GIA uncertainties. Despite similarities between some of the methods shown here, the interpretation of some signals can still vary across methods for the same GIA model. As discussed above, the use of 3-D Earth models may also increase the sensitivity of the results to the starting model. Finally, GIA modeling studies often adjust the forward GIA model to decrease χ 2 values without explicit discussion of model uncertainty. However, studies should at least endeavor to quantify GIA model uncertainty since inferred σ m values are often as large as, or larger than, σ d values.

Summary Points and Recommendations
1. GIA model uncertainty is complex and may take the form of within-model uncertainty, across-model uncertainty, and methodological uncertainty (section 3); modelers and users should be aware that all may affect the ability to separate non-GIA from GIA signals and the interpretation of sea level and mass change budgets. 2. The ±20% rule may underestimate uncertainties in deglaciation centers (and be inappropriate for application in farther-field regions and regional studies). 3. Until improved GIA models are available, a pragmatic approach is to assess what GIA model uncertainties are appropriate for a chosen model and region (e.g., as in sections 4 and 5), provided that adequate a priori non-GIA information is available; the implementation of such an approach precludes the standard χ 2 calculation in which the goodness of fit of the model is evaluated a posteriori since the model uncertainties have been estimated to provide a good fit (χ 2 = 1).