We note that the optimal fits for each region were obtained with different NAIC histories. While this is unrealistic in the sense that there can be only one “true” ice history, we permit this possibility given that we do not know this history and thus apply a suite of models to sample the uncertainty in this model input. The result that different ice models are preferred by data from different regions is due to differences in the sensitivity of the RSL data in each region to aspects of the ice history. For example, data in Canada are sensitive to detailed characteristics of the loading history in that region, whereas data from the Gulf coast are more sensitive to the total volume of ice and its general distribution across North America. Our aim here is to identify model parameters that provide a best‐estimate of the GIA signal and bound this signal to a given degree of confidence. This aim is best served by allowing the possibility of different ice and earth models given our imperfect knowledge of these model inputs.

To partly compensate for the above‐described limitations of the 1D Earth model, we estimated optimal viscosity parameters for three geographic subsets: (i) the two Canadian sites close to the ice margin (Newfoundland and Nova Scotia), (ii) the US Atlantic coast sites, and (iii) US Gulf Coast sites. Other subsets were considered (e.g., including some sites in the NE US with the two Canadian localities) but the results were less optimal. Sub‐division into more than three regions was not considered. Results for these data subsets are shown in Section 3 of Appendix S1 (Figures S10–S12). The results indicate that sensitivity to variations in the ice history diminishes as distance from the location of the NAIC increases. From these results, an optimal model of LT = 120 km, UMV= 2 × 10 21 Pas, LMV = 90 × 10 21 Pas, NAIC model 8034 (ISMI = 27) can be identified for the Canadian data; LT = 71 km, UMV = 3 × 10 21 Pas, LMV = 30 × 10 21 Pas, NAIC model 9894 (ISMI = 34) for the US Atlantic coast dataset; and LT = 120 km, UMV = 5 × 10 21 Pas, LMV = 50 × 10 21 Pas, NAIC model 1243 (ISMI = 2) for the US Gulf Coast dataset.

The results based on the locally scaled model are perhaps the most useful, as they indicate the effect of lateral structure when scaled relative to the center of the study region and not to a global average. These results indicate that the influence of lateral structure on RSL is significant, particularly in the early to mid‐Holocene. Since the majority of the observations are mainly mid‐to‐late Holocene, where the lateral viscosity perturbations result in only a few meters of RSL change, the 1D model is still able to produce quality fits at most sites. What is perhaps more relevant here is the influence of lateral structure on present‐day rates, which is a better measure with regard to the accuracy of the 1D model to predict the GIA contribution to future RSL. Using the results in Figure 3 , the mean of the difference (3D locally scaled model output minus 1D model output) is −0.1 mm/a and the standard deviation is 0.3 mm/a, indicating that the overall bias is small but there is significant inter‐site variability which limits the accuracy of a 1D model. As indicated in Figures 3 and S4, the greatest differences occur in locations proximal to the ice margin. These results are compatible with previous studies that considered the influence of lateral structure on rates of sea‐level change at tide gauge stations [ Kendall et al., 2006 ; Davis et al., 2008 ]. Our 3D model results also provide a partial explanation of why analyses based on 1D models (including this study) have not been capable of producing high quality fits at all locations along the U.S. Atlantic coast [ Engelhart et al., 2011 ; Roy and Peltier, 2015 ]. It is important to re‐iterate that our results are based on only one realization of lateral structure at a relatively low resolution. Further work is required that considers uncertainty in lateral structure to test the accuracy of the preliminary results presented here.

Figures S1–S3 indicate the lateral variation in viscosity structure relative to the optimal 1D profile in the region of interest. These figures show that there are substantial variations in viscosity across the study region within the upper mantle and shallow lower mantle. The influence of this lateral structure on RSL is indicated in Figure 3 , which shows that the largest effects occur at the most northerly sites where the isostatic deformation is the greatest and there are large gradients in the lateral structure. At these sites, structure within the upper mantle dominates the response. At more southerly sites, the influence of lateral structure in the upper mantle tends to produce an RSL signal that is of opposite sign to that produced by structure in the lower mantle, resulting in a net signal that is relatively small (few m). The part cancelation of these signals is broadly compatible with the departures in viscosity relative to the 1D background model, with UMV values being shifted higher and LMV values lower on average (see Figures S1–S3).

To explore the influence of lateral variations in viscosity structure, we utilized the 3D Earth model described above. Predictions were generated for the optimum parameter vector: ice model 1259 of Tarasov et al. [ 2012 ] and a background radial viscosity profile identical to that of the optimal model defined above. The lateral variations superimposed on this 1D structure are based on the seismic tomography model S40RTS; a shear wave velocity model with lateral resolution defined by a spherical harmonic truncation of degree and order 40 [ Ritsema et al., 2011 ]. These lateral perturbations were applied globally throughout the entire depth extent of the mantle as well as only in the upper and lower mantle regions to examine the relative importance of lateral structure at different depth extents. In addition, we chose one of our best fit sites near the center of our region of interest and perturbed the model such that at each depth the viscosity value at our chosen reference site matches that of the 1D background model (referred to as the “Locally Scaled” model).

The lower frame in Figure 2 indicates that relatively low δ values can be produced for UMV/LMV values around 0.3/2 ×10 21 Pas. By considering the δ values for geographic subsets of the data (discussed later in this section; results shown in Figures S10–S12), this local minimum is associated with the US Atlantic coast data (Figure S11). Inspection of Figures S11 and S12 indicates that, as well as having a preference for relatively high UMV and LMV viscosity values, data from the US Atlantic and Gulf coasts are also fit well by some lower viscosity values. The existence of more than one minimum in the data‐model misfit measure is relatively common in GIA analyses [e.g., Milne and Peros, 2013 ; Lambeck et al., 2014 ]. For the entire SLIP database, the parameter vector that produces the best data‐model fits is LT = 71 km, UMV = 3 × 10 21 Pas, LMV = 90 × 10 21 Pas, and NAIC model 1259 (ice sheet model index (ISMI) = 5). While this parameter vector is one of many that provides a good fit to the entire data set (as evident on Figure 2 ), we adopt it as the optimal set as it minimizes the misfit criterion (Equation 1 ).

Data‐model misfit plots for the full dataset as a function of (from top to bottom) upper mantle viscosity and ice sheet model index (ISMI) (lower mantle constant at 30 × 10Pas), lower mantle viscosity and ISMI (upper mantle viscosity constant at 3 × 10Pas), and upper mantle viscosity and lower mantle viscosity (ISMI set to 34; 9894 ice history of]). Note that LT is 71 km for all results shown. Equivalent versions of this figure but for the three sub‐regions defined in the main text are in Figures S10, S11, and S12. These plots show cross sections of a multidimensional parameter space and so do not necessarily indicate all the minima and maxima of the misft parameter.

The misfit parameter δ is shown as a function of key model input parameters for the entire SLIP data set in Figure 2 . The lower frame shows δ as a function of LMV and UMV for a specific ice history and LT value. In this case, the ice model is one that provided some of the best data‐model fits out of the 35 considered. These results indicate that the RSL reconstructions show a preference for relatively high UMV and LMV values. This is also evident when considering the fits as a function of ice history and UMV or LMV. The middle frame in Figure 2 shows δ as a function of LMV and ice model for specific values of LT and UMV; the results indicate that good fits are obtained only for LMV greater than about 20 × 10 21 Pas for the chosen values of UMV and LT. In a similar manner, the top frame indicates that good fits are obtained only for UMV values greater than 10 21 Pas for the chosen values of LMV and LT. These viscosity values are higher than those commonly estimated using near‐field [e.g., Mitrovica and Forte, 2004 ] and far‐field data [e.g., Lambeck et al., 2014 ]; however, they are compatible with a recent estimate based on RSL data from the circum‐Caribbean region [ Milne and Peros, 2013 ].

Data‐model fits were determined by computing, for each parameter vector, a RSL curve at each SLIP location (indicated by colored circles and diamonds in Figure 1 ) and then determining the point on the model curve that is closest to the SLIP. That is, model curves were calculated at each of the almost 700 SLIP locations; the site averaged locations (white crosses in Figure 1 ) are used only for illustrating data‐model fits once optimal parameters are determined. The quality of fit was quantified by squaring the temporal and vertical misfit normalized by the related data uncertainty []. The total data‐model misfit for a given model run is computed by summing up this quantity for all of the SLIPswhere RSLandare the RSL and age values of theth SLIP. RSLandare the values from a given run of the GIA model which is closest to theth SLIP.andare the 2uncertainty values for the RSL and age of the SLIP, respectively. Note that the upper and lower limiting data, although present in Figure 4 , are not included as part of the calculation ofand play no role in the selection of the best fitting model parameters.

2.3.2 Data‐Model Fits and Estimating Uncertainty Bounds

In Figure 4, we show fits for the optimal parameter set that minimized the data‐model misfit (Equation 1) as well as estimates of the 1−σ confidence bounds for each of the three geographic subsets. The limiting data were not used in determining the optimal or bounding model parameter vectors. Estimating the uncertainty range is non‐trivial. Bayes' theorem offers a rigorous approach to uncertainty estimation. However, Bayesian approaches have several requirements in order to attain accurate estimation. These include a well‐defined prior distribution, an accurate likelihood function that fully accounts for structural uncertainties in the model, and a large sampling that covers “reality.” These requirements are rarely fully met for complex physical systems such as that simulated here. We judge that these requirements are too broken in our case and have therefore implemented a more heuristic approach (though still probabilistically self‐consistent) which is described in the following. For the sake of comparison, we also estimate the uncertainty range using a nominally Bayesian approach (“nominally” in the sense that key criteria are not adequately met); the details of which are provided at the end of this subsection.

We estimate one‐sided bounds (one‐sided in the sense that the lower and upper bounds are defined independently) by isolating the model runs that bound, from above or below, (P+(100‐P)/2)% of the SLIP uncertainty boxes, where P is the adopted confidence level (in the SLIPs and the estimated model bounds). For example, (P+(100‐P/2)% is the probability that the actual RSL was, in the case of the lower bound determination, within or above the P uncertainty box of any sample (this assumes that the true RSL is equally likely to be above or below the uncertainty box). We first sought to determine model bounds at the 2 − σ confidence level (P = 95.45%) but no model runs passed in this case. A number of model runs passed at the 1−σ or 68.27% confidence level. Of the subset of runs that passed, the one that provided the best fit to the data was chosen to define the 1−σ bounding runs. To partially address the lack of temporal independence of the SLIPs, the bounding criterion was applied separately to all dates older and younger than 4.5 ka. This method produces physically self‐consistent bounds as they are individual model runs as opposed to an ensemble statistic. It thereby accounts, in part, for the multivariate correlative structure of the misfits that is often ignored in approaches that are otherwise more statistically rigorous due to the difficulty in their ascertainment. The bounds account for uncertainty in the observations as well as parametric uncertainty in the input viscosity and ice history models; however, they do not account for uncertainty associated with model structural error (e.g., assuming 1D viscosity structure). The bounds are conservative in that larger sampling of model parameters (especially over ice chronologies) would likely reduce their range. Computational use of more accurate error ellipses [Parnell and Gehrels, 2015] as opposed to error boxes would likely provide somewhat tighter (and more accurate) bounds but are non‐trivial to implement. Using the above‐described method, two sets of model bounds were estimated: one using all the SLIPs from each region (black dash‐dotted line in Figure 4) and one using only the regional compaction free data subsets (blue dash‐dotted line). The former is associated with the optimal model (red line) since this was determined using all SLIPs from each region. The model parameters associated with the optimal and bounding runs (heuristic only) are given in Table S1.

From Figures 4 and S9, it is evident that the fit for each location is generally of high quality, with the optimal model RSL intersecting the height‐time distribution of SLIPs relatively well. There are a few locations where the optimal model lies below much of the data (e.g., Southern North Carolina) and others where the curve is higher than the majority of the SLIPs (e.g., East Maine). This reflects the fact that the model cannot capture fully the spatial variability in the observations due to model limitations (e.g., spherically symmetric Earth model) or processes other than GIA being significant in some locations (e.g., tectonics).

Figure 4 Open in figure viewer PowerPoint σ time and height uncertainties. We show these rather than the more conventional 2 − σ ranges because no bounding models were found using this more conservative range (see main text for details) and so the bounding runs (dash‐dotted lines) represent 1−σ uncertainty estimates. The black dash‐dotted lines define the estimated 1−σ bounds using the heuristic approach (see Section σ bounds when using the same approach but only compaction free SLIPs in each sub‐region are considered. For comparison, the green dash‐dotted line shows the 1−σ bounds estimated using the nominally Bayesian approach (see Section Sea‐level index points (SLIPs) for selected sites along both the Atlantic coast and the Gulf Coast of North America, along with the optimal (best‐fitting) model time series (red line) for selected locations in each of the three sub‐regions defined in the main text. The SLIPs drawn in blue can be considered compaction free (indicated by “CF SLIP” in key; see Section 3 of Appendix S1). The dimensions of each SLIP box indicate the associated 1−time and height uncertainties. We show these rather than the more conventional 2 −ranges because no bounding models were found using this more conservative range (see main text for details) and so the bounding runs (dash‐dotted lines) represent 1−uncertainty estimates. The black dash‐dotted lines define the estimated 1−bounds using the heuristic approach (see Section 2.3.2 ) when all SLIPs in each sub‐region are considered. The blue dash‐dotted lines define the estimated 1−bounds when using the same approach but only compaction free SLIPs in each sub‐region are considered. For comparison, the green dash‐dotted line shows the 1−bounds estimated using the nominally Bayesian approach (see Section 2.3.2 ) and compaction free SLIPs (indicated by “NB‐CF Bounding” in key). The optimal run (red line) was determined using all the SLIPs in each sub‐region. Limiting data are indicated by either green “T” symbols (upper limiting) or blue inverted “T” symbols (lower limiting). These were not used in estimating the optimal or bounding models but serve as a useful check on the accuracy of these model estimates (see text for further discussion). Note that all model curves (optimal and bounds) are generated at the locations indicated by white crosses in Figure 1 . Results for all 28 of these locations are in Appendix S1 (Fig. S9).

The optimal model for each of the three regions (red lines in Figures 4 and S9) was determined using all SLIPs from each region. To investigate whether downward displacement of SLIPs due to sediment compaction [Törnqvist et al., 2008; Horton and Shennan, 2009] may have influenced the parameter estimates, we also determined optimal parameter vectors using only the compaction free data from the Atlantic and Gulf coasts. This was not done for the Canadian region due to the limited number (8) of compaction free SLIPs. In general, compaction was found to have little effect on our estimates of optimal parameters for each region, resulting in the GIA component of the RSL projections changing by only a few percent. Thus all of the optimal model curves indicated in Figures 4 and S9 are based on all of the SLIPs from each region. However, the estimated model bounds—particularly the lower bound—were significantly affected when considering the compaction free data subsets (compare dash‐dotted blue and black lines in Figures 4 and S9). The greater sensitivity of the bounding runs to changing the input data in this way reflects the conservative nature of the criterion used in determining these runs.

While the limiting data were not used in estimating the optimal and bounding runs, they provide a useful check to determine which of the estimated bounds are most accurate. The lower bounding run for past sea‐level is arguably the most important as it determines the highest RSL rate and thus, the value of the upper uncertainty bound for future RSL change. It is evident from Figures 4 and S9 that the lower bounding run based on all SLIPs violates many of the lower limiting data in each region. For this reason, we believe that the bounds estimated using the compaction free SLIPs are more accurate. The lower bounding runs based on the compaction free datasets are also incompatible with a number of the lower limiting data and so we considered a test to determine whether this incompatibility holds in a statistical sense. Specifically, we took the subset of model runs that passed the SLIP test described above (for compaction free data only) and then applied the following criteria: no more than ((100 − P)/2)% of the lower limiting data P uncertainty error boxes, which indicate uncertainty on the height and time measurement only, should lie above a model RSL curve and at least the same fraction should lie below the curve (such that no more than P% can be intersected by the curve). This criterion follows directly from the assumed symmetric uncertainty for the marine limiting sample elevation and age (the latter is a calculational approximation since the age calibration does not make this assumption). The key difference from the SLIP bounding criteria is that in this case the uncertainty box only defines a one‐way bound for RSL as opposed to defining a time‐elevation region that has probability P of being intersected by the true RSL curve. None of the model subsets that passed the SLIP criterion for each of the three regions satisfied this limiting data criterion, pointing to some combination of model ensemble limitations, inconsistencies between the SLIPs and limiting data, and too narrow observational uncertainties.

The upper bounding runs are less sensitive to whether all SLIPs or those that are compaction free only were considered. Examination of the upper bounding RSL curves shown in Figures 4 and S9 shows that, in some locations, the bounding run based on all SLIPs provides the better result in that it bounds more of the SLIPs and is compatible with the limiting data (e.g., New Jersey); whereas in other locations this bounding run is incompatible with a number of the limiting data (e.g., Connecticut). In considering these differences between locations, it is important to emphasize that the bounding runs are determined using the regional data sets and so the quality of the bounds do vary from location to location within a given region (in the same way as the optimal model fit varies). Given the more mixed results in this case and the lower sensitivity of the upper bound to which data set was used (full versus compaction free), we do not provide a preference for which upper bound to use in each region. Model projections are provided in the following section for both sets of bounds so that the end user can decide which is better to use for their location of interest based on the results shown in Figures 4 and S9. We note that, for the Canadian dataset, the upper bound for the compaction free data sits above that for all SLIPs (the opposite is true for most other localities). This likely reflects the limited number of SLIPs considered for this region.

As previously mentioned, uncertainty ranges were also computed using a nominally Bayesian approach. To construct 1−σ bounds we treat our δ as a log‐likelihood to construct relative probabilities. To generate posterior probabilities, we assume a uniform prior (which is invalid for our case) and then implement Bayes' theorem by normalizing posterior estimates so that they all sum to one over the model ensemble. Thus, we can extract a sub‐ensemble of the highest probability parameter vectors that collectively exceed a chosen posterior probability threshold and we take our bounds as the maximum/minimum RSL value at each location and time from this sub‐ensemble. Compared with the heuristic procedure described above, this method bears more similarity to those employed in previous studies and in this study produces narrower bounds more consistent with these studies [e.g., Kopp et al., 2014]. However, it does not meet the requirements detailed above for accurate posterior estimation. In detail, our δ does not capture the multivariate structure of the residuals, nor does it account for structural uncertainty. The prior is difficult to define given that the ice sheet chronologies were from a glacial model calibration that used a single earth viscosity structure. Furthermore, even though we sample over 12,000 parameter vectors, it is unlikely that our sample space covers “reality” given the dimensionality of possible deglacial histories. A more rigorous Bayesian posterior estimation is possible but well beyond the scope of this study. This would entail some simulation to model the multivariate structure of the residuals, analysis with a greater number of 3D Earth models and glacial systems models to parametrize the structural uncertainty, and the use of statistical emulators to substantially increase the coverage of the sample [e.g., Hauser et al., 2012]. In addition, the prior distribution of ice sheet chronologies would need to be regenerated from a model calibration that sampled over a plausible range of earth viscosity structures.

The impact of these broken requirements is evident in Figures 4 and S9 where the nominally Bayesian approach, when applied to the compaction free data, generally underfits the SLIP data by comparison. Given the societal and environmental context of our study, some over‐estimation of uncertainty is more appropriate than a similar amount of under‐estimation. We therefore focus on the more conservative heuristic approach. While heuristic, the approach is consistent with the stated probabilistic assumptions. It identifies actual model runs (i.e., probabilistic samples) that are consistent with 1−σ confidence intervals. This is evident in the RSL plots if one keeps in mind the requirement that a single run defines a bound for each large region (Atlantic Canada, US Atlantic coast, US Gulf Coast).