Earth and Planetary Science Letters Constraints on mantle viscosity and Laurentide ice sheet evolution from pluvial paleolake shorelines in the western United States

The deformation pattern of the paleoshorelines of extinct Lake Bonneville were among the ﬁrst features to indicate that Earth’s interior responds viscoelastically to changes in surface loads (Gilbert, 1885). Here we revisit and extend this classic study of isostatic rebound with updated lake chronologies for Lake Bonneville and Lake Lahontan as well as revised elevation datasets of shoreline features. The ﬁrst order domal pattern in the shoreline elevations can be explained by rebound associated with the removal of the lake load. We employ an iterative scheme to calculate the viscoelastic lake rebound, which accounts for the deformation of the solid Earth and gravity ﬁeld, to calculate a lake load that is consistent with the load-deformed paleotopography. We ﬁnd that the domal deformation requires a regional Earth structure that exhibits a thin elastic thickness of the lithosphere (15–25 km) and low sublithospheric Maxwell viscosity ( ∼ 10 19 Pas). After correcting for rebound due to the lake load, shoreline feature elevations reveal a statistically signiﬁcant northward dipping trend. We attribute this trend to continent- scale deformation caused by the ice peripheral bulge of the Laurentide ice sheet, and take advantage of the position of these lakes on the distal ﬂank of the peripheral bulge to provide new insights on mantle viscosity and Laurentide ice sheet reconstructions. We perform ice loading calculations to quantify the deformation of the solid Earth, gravity ﬁeld, and rotation axis that is caused by the growth and demise of the Laurentide ice sheet. We test three different ice reconstructions paired with a suite of viscosity proﬁles and conﬁrm that the revealed trend can be explained by deformation associated with the Laurentide ice sheet when low viscosities below the asthenosphere are adopted. We obtain best ﬁts to shoreline data using ice models that do not have the majority of ice in the eastern sectors of the Laurentide ice sheet, with the caveat that this result can be affected by lateral variations in viscosity. We show that pluvial lakes in the western United States can place valuable constraints on the Laurentide ice sheet, the shape of its peripheral bulge, and the underlying mantle viscosity.


Introduction
The western U.S. experienced a mean increase in precipitation during the last glacial cycle, which led to the formation of a series of pluvial lakes that filled the Basin and Range Province (e.g., Benson et al., 1990;Mifflin and Wheat, 1971). The most prominent of those is Lake Bonneville (30-10 ka), an extinct pluvial lake that occupied the eastern Great Basin (Fig. 1B, C). At its maximum extent (∼18 ka) (Oviatt, 2015), the lake had a volume of around Fig. 1. Lake level chronology and shoreline elevations for Lake Lahontan and Lake Bonneville. A, B) Reconstruction of lake level curve for Lake Lahontan (black curve from Benson et al., 2013, which is used here; light blue curve from Reheis et al., 2014) and Lake Bonneville (Oviatt, 2015). For Lake Bonneville, the elevation has been adjusted to account for the rebound of the shoreline (Oviatt, 2015). We extended the Provo stage until 14 ka to test different timings for the duration of the Provo shoreline (dashed line indicates original reconstruction by Oviatt (2015); solid line indicates the curve used here). C) Geographic setting of Lake Lahontan and Lake Bonneville. D-F) Reconstructed still water level (SWL) from the Sehoo, Bonneville, and Provo shoreline features, respectively. Reconstructions are based on the original data by Adams et al. (1999); Chen and Maloof (2017); and Currey (1982). Points that have been removed due to other deformation processes are shown as transparent markers with dashed outline (see SM, Section S1). (For interpretation of the colors in the figures, the reader is referred to the web version of this article.) line features demarcating the lake surface formed on the landscape. After the lakes drained, the solid Earth rebounded, pushing up shoreline features that formed on islands within the deepest part of the lake to elevations higher than those at the lake's margin. Due to this differential uplift in response to the lake unloading, shoreline features at the lake's center are today significantly higher than those on its periphery ( Fig. 1D-F). This pattern is most apparent for features of the Bonneville lake stage, where differences in shoreline feature elevations are over 70 m (Fig. 1E). These paleoshorelines have played an instrumental role in the understanding of isostatic rebound on Earth. Gilbert (1885) reported this phenomenon and provided several possible explanations including isostatic adjustment to the lake load and changes in the gravitational equipotential surface due to load redistributions. While the latter are small (Woodward, 1888), this early assessment paved the way for a gravitationally self-consistent rebound theory that is used in ice age sea level calculations today (Whitehouse, 2018).
The amount of deflection of Lake Bonneville shorelines has been used in numerous studies to constrain Earth's local viscoelastic properties (Bills and May, 1987;Cathles, 1975;Crittenden, 1963;Iwasaki and Matsu'ura, 1982;Nakiboglu and Lambeck, 1982;Passey, 1981). Most recently, Bills et al. (1994) performed an inversion for a multilayer viscosity model and a fixed lake load history resulting in a viscosity profile that exhibits a very low viscosity channel (4 × 10 17 Pa s) beneath the lithosphere that increases to 2.5 × 10 20 Pa s at a depth of 150 km. The lithosphere consists of a thin high viscosity layer (2 × 10 24 Pa s from 0 to 10 km depth) followed by an intermediate layer (∼5 × 10 20 Pa s from 10 to 40 km depth). The rebound pattern of Lake Lahontan was recognized significantly later than that of Lake Bonneville (Mifflin and Wheat, 1971), likely due to its smaller magnitude, with differences in shoreline feature elevations of only up to 25 m (Fig. 1D). Adams et al. (1999) investigated the shoreline feature elevations, and Bills et al. (2007) used them together with a fixed lake load history  (Peltier et al., 2015). The thick gray line indicates the outline of the ice sheet at 18 ka. Lake Lahontan (west) and Lake Bonneville (east) are shown with a black outline, positioned roughly at 40 • N. B) Schematic illustration of the effect of the Laurentide ice sheet on paleolakes in the western United States. Paleoshorelines of lakes on the distal side of the peripheral bulge are predicted to dip down towards the ice sheet today. This is a result of the combined effects of solid Earth deformation, which leads to a downward dip, and a changing gravitational pull of the Laurentide ice sheet, which acts to reduce the total downward dip.
to identify a very low sublithospheric viscosity (less than 10 18 Pa s between 80-160 km). This minimum viscosity is comparable to the values obtained for the Lake Bonneville region but might extend over a larger depth range (Bills et al., 2007). Post-seismic studies from Lake Lahontan find slightly lower viscosities over the same depth range but overlap within uncertainty with the lake rebound obtained viscosities (Dickinson et al., 2016).
These Basin and Range sub-lithospheric viscosity estimates are significantly lower than global average estimates at this depth of ∼5 × 10 20 Pa s, obtained from observations derived from postglacial rebound (Peltier et al., 2015). Both Bonneville and Lahontan lie within the Basin and Range Province, which formed as a result of extension-related faulting (Sonder and Jones, 1999). Joint inversions of seismic and petrologic studies indicate that this region is characterized by a thin crust (30-35 km), shallow lithosphereasthenosphere boundary (50-55 km), and a high asthenospheric potential temperature of 1525 • C (Lekić and Fischer, 2014;Plank and Forsyth, 2016). These elevated sublithospheric temperatures are consistent with body wave tomography results that reveal relatively high S-and P-wave speeds and a high P to S-wave speed ratio, which suggests the presence of sublithospheric melt (Schmandt and Humphreys, 2010). The low viscosity estimates derived from lake rebound studies is therefore consistent with the notion that Earth structure underneath the western U.S. is significantly weaker than cratonic sites such as the Canadian and Fennoscandian shields from which rebound-based estimates of viscosity are normally obtained .
Even after the lake rebound signal is corrected for, longer wavelength spatial trends in shoreline elevations remain. For Lake Bonneville, this residual has largely been attributed to tectonic and crustal deformation such as displacement along the tectonically active Wasatch fault, which straddles the eastern flank of the paleolake (Bills et al., 1994;Nakiboglu and Lambeck, 1982). Similarly, a northward dipping trend in the residual Lake Lahontan shorelines has been linked to tectonics associated with the Yellowstone hot spot (Bills et al., 2007). An alternative explanation put forth earlier by Bills and May (1987) explained a possible northward dipping trend in the residual shoreline of Lake Bonneville with the lake's location on the peripheral bulge of the Laurentide ice sheet. Postglacial rebound calculations of the North American pe-ripheral bulge place these western U.S. lakes on the ice-distal side of the bulge. This long wavelength trend in topography is sampled by these much smaller lakes ( Fig. 2A), resulting in paleoshoreline features that are expected to dip downward towards the ice sheet (Fig. 2B). In addition to solid Earth deformation, the Laurentide ice sheet also deforms the gravity field, exerting a gravitational pull on water in the lake. This effect by itself would cause an upward dip (towards the ice sheet) in paleoshoreline features, counteracting to some extent the downward dipping signal associated with the solid Earth (Fig. 2B).
In this study, we revisit this classic rebound problem to investigate the putative northward dipping trend in the paleolake shoreline features of Lake Lahontan and Lake Bonneville. We use revised shoreline feature elevations and updated chronologies of lake level histories together with state-of-the-art isostatic adjustment modeling to test the hypothesis that a statistically significant northward dipping trend can be detected in all lake stages once the lake rebound pattern is corrected for. We further use three different ice sheet reconstructions together with an ice age sea level calculation to investigate which ice sheet-mantle viscosity structure combination best reproduces the observed lake tilt. While proglacial lakes have been used to constrain ice sheet evolution and mantle viscosity in recent work (Gowan et al., 2016;Lambeck et al., 2017), this study is the first to investigate the deformation of distant pluvial lake shorelines.

Lake chronology
Lake Bonneville and Lake Lahontan were the two largest pluvial lakes in the Great Basin during the last Pleistocene glaciation (∼30-10 ka; Fig. 1C). A common misconception is that these lakes were hydrographically connected to the Laurentide or Cordilleran ice sheets as ice-dammed or glacial lakes; in actuality, these lakes were hydrographically distinct and instead fed by local precipitation and snowmelt delivered by perennial rivers (Reheis et al., 2014). Both lakes occupied basins of similar topographic characteristics, filling in broad and flat valley floors surrounded by steep mountainsides consistent with the extensional tectonic regime of the Basin and Range Province. Although both lakes were essentially contemporaneous, the lake level history of Lake Bonneville is better constrained. Lake Bonneville was a deeper lake existing as a single entity over a greater period of its history. Thus, a reconstruction of its lake level history that is consistent with most interpretations of sediment core and outcrop evidence has been more feasible. In contrast, the shallower Lake Lahontan existed as several smaller, disconnected sub-basins for most of its history, complicating attempts to reach consensus on its lake level history (Benson et al., 2013;Bills et al., 2007;Reheis et al., 2014).
Figs. 1A and 1B depict the most recent reconstructions of lake level histories for Lake Lahontan and Lake Bonneville, respectively. Both histories were derived from many decades of extensive field and sediment core observations and are constrained by hundreds of radiometric dates on organic material, tufas, and tephras extracted from cores, exposed outcrops of lake deposits, and deposits associated with geomorphic shoreline features (e.g., Adams et al., 1999;Briggs et al., 2005;Oviatt, 1997;Oviatt et al., 1992;Patrickson et al., 2010;Sack, 2015;Spencer et al., 2015). During their existence, each lake experienced a rise (transgression) and fall (regression) of water level, and in certain instances, left behind evidence of their evolution in the form of prominent shorelines features. Of the many sequences of paleoshorelines available, the shoreline features most relevant to this study are those associated with the maximum extent of Lake Lahontan, the Sehoo lake stage (∼15 ka; Figs. 1A and D), and the Bonneville and Provo lake stages of Lake Bonneville (18 ka and 18-15 ka; Figs. 1B, E and F). Evidence suggests that Lake Bonneville did not occupy its maximum extent, the Bonneville lake stage, for more than a few hundred years (Gilbert, 1885;Oviatt and Jewell, 2016) before a catastrophic collapse of an alluvial-fan dam dropped lake levels by 100 m to settle at the Provo level (Miller et al., 2013).
For clarity and simplicity, we hereafter use Sehoo in reference to the stage at which Lake Lahontan reached its greatest extent (e.g., the Sehoo shoreline or Sehoo lake stage), and Bonneville and Provo in reference to those stages associated with Lake Bonneville's history (e.g., the Bonneville shoreline or Provo lake stage). The phrases Lake Lahontan and Lake Bonneville will only be used when referring to the entire lake cycle, encompassing all fluctuations depicted in Figs. 1A, B; earlier major lake cycles in these basins exist and have other names (Oviatt et al., 1999).
We note that there are differences in the degree of certainty in the different lake level reconstructions. For example, it is thought that the timing of the Bonneville lake stage is much better constrained than the end of the overflowing phase at the Provo shoreline (Oviatt, 2015). In the case of Lake Lahontan, different interpretations of sediment cores and outcrops have also led to conflicting lake level reconstructions (Reheis et al., 2014). Despite these nuances, our experimental design requires that we take the interpreted lake level curves at face value. We use the lake level histories by Oviatt (2015) and Benson et al. (2013) for Lake Bonneville and Lake Lahontan, respectively, to constrain the temporal evolution of the lake load in our model, one of the key initializing inputs in our workflow (Fig. S1). In each iteration, we update the lake level curve such that it coincides with the shoreline feature elevations on the modeled paleotopography (see Section 3.1 and Fig. S1). Lastly, we test the sensitivity of our results to the timing of the end of the Provo lake stage.

Shoreline data
We use elevation data of shoreline features from three sources: Adams et al. (1999), which provides data for the Sehoo shoreline of Lake Lahontan; Currey (1982) for both the Bonneville and Provo stages of Lake Bonneville; and Chen and Maloof (2017) for the Bonneville stage of Lake Bonneville. We note that an important part of the study carried out by Chen and Maloof (2017) was a revisitation of the Bonneville shoreline feature data collected by Currey (1982). Because Currey (1982) carried his study out prior to GPS availability, approximately half of his sites were remeasured with modern differential GPS technology (Chen and Maloof, 2017). Therefore, while we use the Currey (1982) dataset of Provo shoreline features in its entirety, we combine both datasets by Currey (1982) and Chen and Maloof (2017) for our analysis of Bonneville shoreline features, opting to use revisited measurements by Chen and Maloof (2017) when available. In total, these datasets provide shoreline feature elevation constraints at 170 sites for the Sehoo lake stage; 274 sites for the Bonneville lake stage; and 112 sites for the Provo lake stage (Fig. 1D-F).
In order to use all three datasets simultaneously, additional processing is required. First, the longitude, latitude, and elevation data are converted to use the same coordinate system and vertical datum: the North American Datum of 1983 (NAD 83) and the North American Vertical Datum of 1988 (NAVD 88) (see Supplementary Material, SM, for details). Second, we address potential biases introduced by differences in the tools and methods used to measure shoreline feature elevations. While all the data by Adams et al. (1999) and Chen and Maloof (2017) were in-field measurements made by total station survey and GPS, the data by Currey (1982) have been shown to generally overestimate the elevation of features (by 1.8 ± 1.4 m, on average) (Chen and Maloof, 2017). Therefore, we apply an adjustment to the data by Currey (1982) based on the method used for each site (see SM for details). Third, we address potential biases introduced by different shoreline feature types in each dataset. Along a shoreline of a lake, many processes associated with the same body of water can form adjacent shoreline features with differing morphological characteristics. Such features include spits, barrier ridges, pocket barriers, wave-cut terraces, and incised alluvial fans (e.g., Adams and Wesnousky, 1998;Chen and Maloof, 2017). Because we are interested in the solid Earth deformation pattern captured by shoreline features that record the position of the mean formative water surface (the still water level; SWL), we must consider differences in how this surface is manifested by each type of shoreline feature. To account for such differences we implement a scheme similar to that of Chen and Maloof (2017) to determine SWL constraints from elevational measurements of shoreline features gathered by Currey (1982) and Adams et al. (1999) (see SM for details).
Because we are solely interested in understanding the deformation pattern induced by lake rebound and a possible regional tilt, we also remove from our analysis shoreline features which are known to have, or are strongly suspected of having, undergone a non-negligible amount of local, post-depositional displacement by other processes. Examples of such excluded data include shoreline features associated with the Wasatch Fault flanking the eastern boundary of Lake Bonneville, and localities associated with Pahvant Butte or Cove Creek Dome that have undergone volcanic deformation since the Holocene (see SM for details, Fig. 1E).

Viscoelastic model
We calculate the deformational and gravitational response to Pleistocene lake and ice loads globally using a spectral approach with spherically symmetric Maxwell rheology (Peltier, 1974). Previous work employed a half-space geometry and did not account for gravitational effects (Bills et al., 2007(Bills et al., , 1994.

Lake rebound modeling
Calculating the response of the solid Earth to changes in the pluvial lakes requires inputs of Earth's internal viscoelastic structure and the temporal evolution of the lake load. We perform a  (Peltier et al., 2015), the ANU ice model (Lambeck et al., 2017) and NAICE (Gowan et al., 2016), respectively. D) Sea level equivalent ice volume during the deglaciation for the southwestern part of the Laurentide ice sheet (red box in panel C).
suite of calculations in which we vary the elastic thickness of the lithosphere and sub-lithospheric viscosity. It is important to note that the elastic thickness of the lithosphere, as utilized here, is a quantity that can differ from lithospheric thickness estimates obtained from seismology or geochemistry (Watts et al., 2013). In our calculations, the lithosphere is treated as a completely elastic solid, while the underlying mantle is treated viscoelastically.
The lake volume could be estimated using the elevation of the lake shoreline (Fig. 1A, B) and the present-day topography. However, this approach underestimates the lake volume because it neglects the downward deflection of the lake basin when the lake load was present. For example, for Lake Bonneville, the lake volume would be underestimated by nearly 20% (Fig. S2C). Thus, estimates of paleolake volumes and lake level curves are dependent on the spatial pattern and magnitude of lake rebound, and vice versa. To avoid this circularity, we iteratively calculate the lake volume, self-consistently accounting for the deflection of the solid Earth and its gravity field (Fig. S1).
We begin with an initial estimate of the lake volume that we derive by filling the present-day topography following a given lake level curve (Figs. 1A, B). We use present-day topography from ETOPO1 (Amante and Eakins, 2009) (approx. 1.8 km spatial resolution). Next, we step through time, calculating the gravitational and deformational response to this changing load for a given viscoelastic Earth structure (Peltier, 1974). We do not assume isostatic equilibrium at each timestep but account for the full time-dependent viscoelastic and gravitational response. Since this response is smooth, this calculation can be performed at a coarser resolution (ca. 20 km). The resulting time-varying topographic change is linearly interpolated onto a grid of higher resolution and combined with present-day topography to obtain a time-dependent, reconstructed, high-resolution paleotopography.
The adjusted topography, together with the lake level curve, is then used to re-calculate the time-dependent lake volume. This new lake volume is once more used to calculate the solid Earth response. We iterate over this procedure until the solid Earth response and the lake volume remain unchanged for any further iteration. In each iteration, we aim to verify that the prescribed lake level curve fills the lake up to the observed SWL (and not higher or lower) during the Sehoo, Bonneville, and Provo lake stages. To accomplish this goal, we include one additional step in each iteration. After the deflection due to lake loading is calculated, we determine the adjusted elevation of the observed SWL on the new paleotopography. For example, if the SWL is inferred to be at 1550 m at a certain location today and loading deflected this site down to 1530 m, the adjusted elevation of SWL corrected for deformation is 1530 m. We next update the lake level curve to fill the lake up to the mean adjusted elevation of all SWL data points during the lake stages for which we have observations (Sehoo, Bonneville, and Provo). This iterative procedure results in the three self-consistently calculated quantities: (1) reconstructed paleotopography, (2) lake level histories and (3) lake volumes for both Lakes Bonneville and Lahontan (Fig. S2).

Glacial isostatic adjustment modeling
In order to calculate the response of the solid Earth to the changing Laurentide ice sheet, we use a gravitationally selfconsistent approach to solve the sea level equation (Kendall et al., 2005). This approach takes the redistribution of water between ice and oceans into account and accurately captures the migration of coastal shorelines and changes in Earth's rotation axis. We use three ice models for our ice load: ICE-6G (Peltier et al., 2015), the LW-6 ANU ice model (Lambeck et al., 2017) and NAICE (Gowan et al., 2016). For ICE-6G, we remove mountain glaciers in the western U.S. because their size in this ice model is significantly larger than actual reconstructions of glacier sizes derived from moraine studies (Fig. S3, Table S1 and references therein).
All ice models are based on different sets of relative sea level curves from around the ice sheet and GPS observations of presentday solid Earth deformation. The ANU and NAICE models further use proglacial lake levels as constraints in their reconstruction. The ICE-6G and NAICE models use a fixed Earth viscosity model and invert for the ice evolution, while the ANU model jointly inverts for ice and Earth parameters. The ICE-6G and ANU models consider the requirement of matching the LGM sea level lowstand and adjust the ice evolution without explicit ice physics requirements. In contrast, the NAICE model does not attempt to match the large LGM ice mass needed to match the LGM sea level lowstand, but does employ a physical ice sheet model to determine the shape of the ice sheet.
As a result of the wide variety of constraints adopted, the ice models vary significantly. Notably, the ICE-6G ice sheet has the largest volume overall and the ICE-6G and ANU ice models have large ice domes over Hudson Bay, which is absent in the NAICE model (Fig. 3). The ice volume in the region just north of the western U.S. (red square in Fig. 3C) is also largest in the ICE-6G model and similar in size between the ANU and NAICE model (Fig. 3D). Furthermore, their temporal evolutions differ. Both the ICE-6G and ANU ice models reach their maximum extent early (26 ka), which is related to fitting the LGM sea level lowstand, while the NAICE model exhibits a later maximum ice extent around 19 ka, during which ice mass in the southwestern Laurentide exceeds that of the ANU ice model. The main ice retreat in the NAICE model occurs after 17 ka, while it is later in the ANU and ICE-6G model (after 14.5 ka).
We pair each ice model with different models of Earth's internal viscoelastic structure and compare the resulting shape of the peripheral bulge to the lake rebound corrected shoreline elevations. The formation and collapse of the Laurentide ice sheet's peripheral bulge itself also affects the spatial distribution of the lake load described in Section 3.1. Therefore, we must perform an additional suite of iterative lake rebound calculations in which we include this ice-related deflection in the time-dependent paleotopography that is used to self-consistently calculate the lake volume (Section 3.1, Fig. S1).

Deflection due to lake rebound
We investigate the fit between the reconstructed SWL and the predicted lake rebound that is obtained using the algorithm outlined in Section 3.1. For the Earth model, we vary the elastic thickness of the lithosphere from 10 km to 30 km and explore sub-lithospheric viscosities spanning a range of 5 ×10 18 to 5 ×10 19 Pa s, guided by earlier inversions (Bills et al., 2007(Bills et al., , 1994. The lake rebound pattern will only be sensitive to shallow mantle structure given the limited lateral extent of the lake. We used perturbation theory to calculate the viscosity sensitivity kernel (Lau et al., 2016) and found that for a reference model with a 25 km thick elastic lithosphere and 10 19 Pa s a upper mantle viscosity, the sensitivity of the lake rebound induced surface deflection rapidly decreases below 300 km. We therefore choose the sub-lithospheric viscosity to extend to 300 km depth, and fix the viscosity structure beneath that to the standard VM5 viscosity profile (Peltier et al., 2015). For the elastic and density structure, we assume PREM (Dziewonski and Anderson, 1981). To evaluate the misfit between our predictions and the observations, we use the reduced chi-squared metric, χ 2 red : where n is the number of observations, m is the number of fitted parameters, p i and o i are the i-th predicted and observed SWL, respectively, and σ i is the latter's uncertainty. The smaller this metric, the better the fit, however, once χ 2 red is 1, the fit is as good as would be expected given the uncertainties in the observations. We will report the number of fitted parameters (m) that we use in each calculation throughout the study. Predictions are compared to SWL data at 15 ka for the Sehoo lake stage, 18 ka for the Bonneville lake stage, and 15 ka for the Provo lake stage, unless noted otherwise.
We find that the best fitting Earth model is not the same between the different lakes (Fig. 4). While the Sehoo and Bonneville lake stages are most consistent with an elastic thickness of 20-25 km and a low viscosity (< 2 ×10 19 Pa s), best fits for the Provo lake stage are obtained for a slightly thinner elastic thickness (15-20 km) and a stiffer underlying mantle (> 3 × 10 19 Pa s). The higher misfits at low viscosities for the Provo lake stage are due to an underprediction of the magnitude of rebound. A larger magnitude can be obtained if the lake at the Provo lake stage is not in isostatic equilibrium but instead still experiences remnant defor- Comparison between the predicted SWL and the observed SWL for our preferred Earth model (white box in Fig. 4). Panels A, B, and C show results for the Sehoo, Bonneville, and Provo lake stages, respectively. Underlying contours show the model prediction while overlain circles show the data. Panels D-F show this comparison as a function of latitude. Black markers are observations and their associated uncertainties, red markers are the model prediction. Error bars represent 1-sigma range uncertainties for SWL estimates. mation from the larger magnitude Bonneville deformation. Therefore, the discrepancy between the Bonneville and Provo shorelines can be reduced if the shoreline features of the Provo lake stage formed earlier than 15 ka. Sensitivity tests demonstrate that the region of best fit is pushed towards weaker viscosities for earlier times of formation (Fig. S4). If the Provo shoreline features formed 1,000 years earlier (16 ka), the best-fitting Bonneville and Provo lake stage viscoelastic models would be more consistent. This earlier time of formation is within the data uncertainty of the Provo shoreline (Oviatt, 2015). Another possible explanation for this discrepancy is that a two-parameter Earth model is not sufficient to capture the deformation of the shoreline. Repeating our analysis with the viscosity profile by Bills et al. (1994), which uses a 9layer viscosity profile (including the lithosphere) that was inverted for using similar shoreline elevations, results in χ 2 red values of 0.62, 11.5, and 3.0 (m = 9) for Sehoo, Bonneville, and Provo lake stages, respectively. Note that the χ 2 red for the Bonneville lake stage is higher due to the lower data uncertainty. With exception to the Provo lake stage, these χ 2 red values are higher than those obtained using our best fitting viscosity structures (χ 2 red of 0.61, 10.2, and 3.3 (m = 2) for Sehoo, Bonneville, and Provo lake stages, respectively; viscosity structure marked by white box in Fig. 4).
For the remainder of this study, we will use a model with an elastic lithospheric thickness of 20 km and a sublithospheric viscosity of 2 × 10 19 Pa s (white box in Fig. 4), which reasonably captures the rebounds of Lake Lahontan and Lake Bonneville (Fig. 5). We infer a volume of 10,250 km 3 and 4,920 km 3 for the Bonneville and Provo lake stage, respectively (Fig. S2C), which is in agreement with the Bonneville volume estimates made by Chen and Maloof (2017), but slightly smaller than estimates by Adams and Bills (2016) who obtained 10,420 km 3 and 5,290 km 3 for the Bonneville and Provo lake stage, respectively. The volume for the Sehoo lake stage is 2.0 km 3 (Fig. S2F), which is slightly smaller than the value of 2.2 km 3 used by Bills et al. (2007) that is based on shoreline elevations from Benson et al. (1995).
To investigate any remaining residual deflection in the shoreline data, we remove the predicted lake rebound pattern from the observations (Fig. 6A-C). We find a noticeable north-south trend in the residual shoreline data, tilted down towards north (Fig. 6D-F). We employ the Mann-Kendall test to investigate whether there is a monotonic (upward or downward) trend in the data that significantly differs from zero (Kendall and Gibbons, 1990). For all three lakes, the results of the Mann-Kendall test indicate that there is a north-south trend (99.9% level of significance).

Deflection due to ice peripheral forebulge
Next, we test the hypothesis that the trends detected in the lake rebound corrected shoreline observations (Fig. 6D-F) are caused by the long wavelength deformation associated with the Laurentide ice sheet. We perform a suite of ice age calculations following the method described in Section 3.2. These model predictions will be sensitive to the evolution of the Laurentide ice sheet and viscosities at greater depth compared to the lake rebound given the larger spatial extent of the Laurentide ice sheet. We explore three ice models (ICE-6G, ANU, and NAICE) and a suite of mantle viscosities. To maintain a fit to the lake rebound patterns, we construct viscosity profiles that follow our best fit model from Section 4.2 (20 km elastic lithospheric thickness and 2 × 10 19 Pa s sublithospheric viscosity) and vary the viscosity between 300 km and the base of the transition zone (670 km) from 10 20 Pa s to 10 21 Pa s and the viscosity between the base of the transition zone and 1175 km depth from 5 × 10 20 Pa s to 5 × 10 21 Pa s (gray shaded bands in Fig. 7). Each ice model is associated with a specific viscosity profile (Fig. 7), which mostly represents mantle structure underneath the Canadian shield. We deviate from these profiles here in order to investigate Earth structure underneath the western U.S. and pair each ice model with different models of Earth's internal viscoelastic structure.
Ideally we would like to perform calculations with lateral viscosity variations. However, these calculations are computationally expensive and not well suited to explore the parameter space. We therefore perform calculations with radially symmetric viscosity structures that allow running many ice-viscosity scenarios. However, this approach comes at the expense that viscosity profiles are global and, in this case, incorrect in locations such as Hudson Bay. In light of this, we perform one additional calculation with lateral variations in viscosity to explore this model limitation (see SM).
Once more, we determine χ 2 red between the observations and predictions, which now includes both the prediction for lake rebound and ice peripheral forebulge deformation. In the lake rebound calculation, we now include the ice-related tilt in our paleotopography, causing slight movement of the water load towards the southern part of the lake that modifies the loading. To test whether the fit to the data is significantly improved when a modeled ice-related tilt is included, we use a two-sample F -test. This test assesses the degree to which the variance in the lake rebound corrected observations is distinct from the variance in the observations that are corrected for both the lake rebound and the ice-related tilt, accounting for uncertainties in the observations.

Trends in viscosity
Modeling results show that the higher the viscosity (in parts of the upper or lower mantle), the larger the predicted tilt across the forebulge at the location of Lake Bonneville and Lake Lahontan. Increasing viscosity in the parts of the upper and lower mantle that we vary here, leads to a higher viscosity contrast across 300 km depth and 670 km depth, both of which results in flow that is more localized at shallow depth, leading to a steeper periph- Fig. 7. Viscosity profiles. Profiles that are associated with the different ice models are shown in color. Note that the E-6 ANU model in purple is the best-fitting Earth model for the ice model LW-6 used here and is provided with an uncertainty. The gray bands indicate the range over which we varied the viscosity. Only certain viscosity profiles are permitted by the tilt in the Bonneville shorelines (see Fig. 8). The black viscosity profile corresponds to one of the best fitting profiles for the western U.S. based on fitting the tilt in the lake rebound corrected paleo shorelines (this viscosity model is outlined by a white box in Fig. 8). eral bulge. The sensitivity to the ice-related tilt is largest for the Bonneville lake stage (Fig. 8B, E, H). At high viscosities, the icerelated tilt that is predicted is larger than the tilt obtained from the lake rebound corrected elevations, leading to an increase in χ 2 red compared to no ice-related tilt correction (purple color, Fig. 8). However, at lower viscosities, the ice-related tilt is flatter and results in a good fit to the observed tilt in the lake rebound corrected shoreline observations (green color, Fig. 8). For the Bonneville lake stage, the F -test reveals that the spread of the residuals is significantly improved when low viscosity Earth models are adopted (Fig. 8B, E, solid line 90% significance level; dashed line 85% significance level). The χ 2 red metric shows that for the Sehoo and Provo lake stages, the fit improves for most viscosity structures when the ice-related tilt is considered (especially Fig. 8A, C, D, F), but the spread of the residuals is not significantly reduced (note how no areas are outlined by a black solid or dashed line).

Trends across ice sheet reconstructions
For the ICE-6G ice model, tilt predictions for the Bonneville lake stage match the lake rebound corrected observations for low viscosities in the parts of the upper and lower mantle that are varied here, with trade-offs between the two (black outline, Fig. 8B). The ANU ice model does not lead to a significant reduction in the variability of residuals (at 90% significance) for the Bonneville lake stage, which indicates that, for our viscosity range, this ice model does not capture the tilt as well as the other ice models. Overall, the χ 2 red values vary less between runs for the ANU model, which suggests that the sensitivity to viscosity variations is lower for this ice model. The NAICE model leads to similar results compared to ICE-6G, despite the significant differences in the ice history (Fig. 8B, H). For the Bonneville lake stage, tilt predictions match the lake rebound corrected observations best for low viscosities in the parts of the upper and lower mantle varied here, with trade-offs between the two (black outline, Fig. 8H). Lastly, this ice model shows most sensitivity to the Sehoo and Provo shorelines because it results in the largest peripheral bulge for high viscosities.
The ICE-6G ice model has significantly more ice volume than the other ice models, which leads to a larger peripheral bulge (Fig. S5A). However, more ice volume also results in stronger gravitational attraction, counteracting the tilt in the paleoshorelines caused by peripheral bulge deformation (Fig. 2B, Fig. S5B). Considering the Bonneville lake stage, ICE-6G results in a large peripheral bulge that is only somewhat compensated by the self-gravitation effect of the ice sheet, causing a significant tilt across the lake (Fig. S5C). The smaller NAICE ice model on the other hand leads to a smaller peripheral bulge, but also less self-gravitation resulting in the preservation of the tilt signal across the lake (Fig. S5G-I). In the ANU ice model, the ice distribution in the western Laurentide ice sheet is significantly smaller than the eastern Laurentide ice sheet (Fig. 3B). As a consequence, the peripheral bulge is centered on South Dakota to the northeast of Lake Bonneville rather than directly north as is the case for ICE-6G and NAICE (Fig. S5D-F). Therefore the ANU ice model leads to slightly less sensitivity to the specific viscosity profile and a worse fit to the clear northsouth trend in the residuals.
Considering the lake stages at 15 ka, the NAICE model leads to the largest peripheral bulge, despite significantly less ice volume than the other ice models (Fig. 3D). This result can be explained by the interplay of forebulge deformation and self-gravity of the ice sheet. In the NAICE ice model, the ice sheet retreats rapidly after 17 ka associated with the collapse of the Laurentide-Cordilleran Ice Sheet saddle. In response to this retreat, the peripheral bulge slowly subsides, leading to a peripheral bulge at 15 ka that is smaller than the one in both ICE-6G and ANU, but still significant. The loss of self-gravitation associated with the ice sheet is, in contrast, instantaneous, and thus no longer counteracts the deformation associated with the peripheral bulge. The combined result is that NAICE exhibits the largest tilt signal among the three ice models (Fig. S6).

Discussion of preferred ice-Earth model
The best fitting ice-Earth model (lowest overall χ 2 red ) is the NAICE ice sheet model paired with a viscosity of 2.5 × 10 20 Pa s between 300-670 km depth and 5 × 10 20 Pa s below that (white rectangle, Fig. 8H). A direct comparison between the lake rebound corrected shorelines and the ice-related tilt from this model shows good agreement (Fig. 9). As described in Section 3.2, this calculation includes a recalculation of the lake rebound that takes the icerelated deformation (solid Earth tilt and gravitational effects) into account, which causes the distribution of the water load to shift southward. This adjustment leads to less rebound in the northern part of the lake and slightly more rebound in the southern part, resulting in deformed contours within the lake (Fig. 9A-C). Overall this process acts to slightly decrease the inferred water volume for Lake Bonneville resulting in a volume of 10,187 km 3 and 4,893 km 3 for the Bonneville and Provo lake stage, respectively. After the correction for the ice-related tilt, there is no longer a significant north-south trend in the corrected observations ( Fig. 9D-F, significance level 95%). Using our iterative approach to calculating the lake rebound and tilt corrected paleotopography, and selfgravitational effects of the ice sheet, we provide gridded datasets of reconstructed water depth for the Sehoo, Bonneville, and Provo lakes (see section 'Datasets', Fig. S7).
Our best fitting Earth structure models have viscosities that are low relative to a Laurentide-centered viscosity model throughout the mantle (Fig. 7). However, trade-offs exist and may allow for higher viscosities in the lower mantle, which would Misfit between the residuals (SWL corrected for lake rebound) and the ice-related tilt calculated from different ice and Earth models. Panels A-C, D-F, G-I show results for the ICE-6G, ANU, and NAICE ice model, respectively. Left, middle, and right panels show results for the Sehoo, Bonneville, and Provo lake stages, respectively. The viscosity structure varies in parts of the upper mantle (between a depth range of 300-670 km, vertical axis) and lower mantle (between a depth range of 670-1175 km, horizontal axis). Above 300 km the Earth structure is identical to what is used in Figs. 5 and 6. The viscosity below 1175 km is 3 × 10 21 Pa s. The misfit is quantified as χ 2 red (Eq. (1) with m = 4) and the color scale is centered on the χ 2 red value obtained without a correction for the ice-related tilt. Purple colors indicate that the fit is worse when the ice-related tilt is accounted for, whereas green colors show that the fit improves. Tiles outlined in black indicate runs that show a significant improvement when the ice-related tilt is corrected for (based on the F -test, solid line is 90% significance level, dashed line shows 85% significance level). The white box indicates the model parameters used in Figs. 9 and 10 and shown by the black line in Fig. 7. require an even lower viscosity in the upper mantle between 300-670 km depth. The low viscosity throughout the upper and lower mantle and the corresponding muted deformation of the peripheral bulge is consistent with sea level indicators along the U.S. West Coast (Creveling et al., 2017). A low viscosity across the upper mantle has also been found for far-field sea level sites (Lambeck et al., 2017) and underneath the Amundsen Sea (Barletta et al., 2018). However, at greater depths (> 400 km), seismic tomography suggests the presence of slab fragments associated with multiple stages of subduction (Sigloch, 2011), which we would expect to result in higher viscosities compared to what is found here. Comparison between the residuals (SWL corrected for lake rebound) and prediction from ice age calculation for the viscosity model shown outlined in white in Fig. 8. Panels A, B, and C show results for the Sehoo, Bonneville, and Provo lake stages, respectively. Underlying contours show the model prediction, while circles show the observations. Note the deflection of contours within the lake that arise from additional lake loading when the ice-related tilt is accounted for in the lake rebound calculation. Panels D-F show the residuals after correction for the ice-related signal as a function of latitude. Marker sizes in all panels are inversely proportional to the data uncertainty. Uncertainties are scaled the same in panels D and F, but are different in panel E (applying the same scaling would lead to very large markers).

Limitations of this analysis
There are several limitations to our results. First, our results are non-unique and trade-offs exist in the viscosity model and the ice sheet reconstruction. It is likely that including additional, potentially high viscosity intermediate layers interspersed with lower viscosity layers could produce an equally good fit to the observed ice-related tilt. Trade-offs also exist in the ice sheet reconstruction regarding the size of the ice sheet, the time of ice growth and the spatial distribution.
Second, while the observations derived from the lake rebound process represent a local constraint on subsurface viscosity structure, the ice-related tilt has sensitivity to viscosity structure that extends to depth (as explored here) as well as laterally, towards the former ice sheet . We explore this sensitivity with one additional exploratory simulation (see SM, Figs. S8 and S9) and find that lateral variations in viscosity can affect the direction of the peripheral bulge tilt. Particularly we find that low viscosities associated with the Yellowstone hotspot can lead to a northeast-southwest tilt in the prediction. While this trend is not evident in the data, when combined with the ANU ice model it could lead to a more north-south trending forebulge. These results are very sensitive to the specific shape and magnitude of lateral viscosity variations, which remain poorly understood. Exploring these, paired with a variety of ice models, requires more efficient inversion schemes for models with lateral variability in Earth structure, which are currently being developed .

Remaining patterns in shoreline elevations
While the ice-related tilt can explain a significant portion of the lake rebound corrected elevations, systematic residuals persist (Fig. 10). The χ 2 red parameter is below 1 for the Sehoo lake stage, which suggests that the shoreline elevations can be explained by our two modeled processes, within observational uncertainty. By contrast, χ 2 red remains above 1 for Lake Bonneville, indicating that additional mechanisms of post depositional deformation are required to explain the spread in the data. Particularly, there is an additional east-west trending pattern in the lake rebound and ice-related tilt corrected shoreline elevations of Lake Bonneville (Fig. 10A, B). The eastern flank of Lake Bonneville is bordered by the Wasatch fault, which has been active since the formation of the paleolake shorelines (USGS, 2017) and vertical displacements since the Holocene are on the order of meters (DuRoss, 2008). Additional parallel faults exist that have experienced less displacement (Fig. 10C) (Friedrich et al., 2003). The pattern of low residuals on the WSW and ENE side, and high residuals in a NNW-SSE strip down the middle could be associated with NNW trending tilted fault blocks or a cylindrical fold associated with continuing tectonic activity on the Wasatch and parallel faults. A comparison of the residuals to the locations of deltaic depocenters (Currey, 1982) (2017), glacial ice caps from Laabs and Munroe (2016), and sediment depocenters from Currey (1982). and glacial ice caps (Laabs and Munroe, 2016) that might have caused additional deformation does not reveal any obvious spatial relationship (Fig. 10C).

Conclusions
We revisit the deformed elevational pattern of Lake Bonneville and Lake Lahontan shoreline features to investigate the different contributions to their deformation. The first order signal is the unloading of these extinct lakes, which leads to a domal deformation pattern in the lake shorelines. In line with previous work, we find that the degree of lake rebound is indicative of a thin elastic lithosphere and weak upper mantle, consistent with the wider tectonic context of this region. The viscosity structure inferred from the Provo and Bonneville lake stage data are most consistent with each other if the Provo shorelines formed early (16 ka). Upon correction for lake rebound, we find that the residual shorelines show a systematic and statistically significant northward dipping trend, a pattern that is consistent with a regional tilt induced by the peripheral bulge created by the extinct Laurentide ice sheet. We perform a suite of ice age calculations and find that the fit to the shoreline data is improved when we include the loading and associated deformation of the Laurentide ice sheet. We explore what ice sheet reconstructions and viscosity profiles produce the best fit to the observed shorelines. We find that while ice volume is a primary control on the size of the peripheral bulge, this effect is counteracted by self-gravity of the ice sheet, resulting in a good fit between the rebound corrected shoreline observations and the predicted tilt for both large (ICE-6G) and small (NAICE) ice sheets. However, the ice distribution affects the size and orientation of the peripheral bulge and we find that an ice model with most of its ice volume in the eastern Laurentide (ANU) is less compatible with the rebound corrected shoreline observations. Lateral variations in Earth's viscoelastic structure can also affect the orientation of the peripheral bulge and might counteract this misfit. Largely independent of the ice sheet model, we find that the tilt is only obtained when the viscosity profile exhibits low viscosities relative to Laurentide centered estimates, which could occur in the upper or lower mantle. Since this result is consistent across the different ice models, it supports the emerging notion that lateral variations in Earth's internal properties are significant and must be considered in global sea level studies (Li et al., 2018). Remaining residuals likely are related to tectonic deformation with possible implications for seismic hazard assessment.

Author contributions
JA led the experimental design of the work, implemented the model simulations, and directed the analysis and interpretation of results, with input from CYC, HL, and ACM. CYC and ACM initiated the project. CYC gathered and reassessed the observational data, with input from JA. KL ran the 3D finite volume model simulations, with input from JA. JA led the drafting of the figures and manuscript text; CYC contributed text regarding the observational data. All authors made significant contributions to the editing of the manuscript.

Datasets
Datasets related to this article can be found at https://doi .org / 10 .5281 /zenodo .3576251 on Zenodo, an open science and open data repository. Datasets include updated shoreline feature elevations for the Bonneville and Provo lake stages of Lake Bonneville; paleotopography (deformed due to lake and ice sheet loading) and lake water-depth reconstructions of Lake Bonneville and Lake Lahontan from 30.5 ka and 29 ka, respectively, to present-day.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.