The influence of glacial isostatic adjustment on continental shelf stratigraphic correlation

Continental shelf stratigraphic architectures are unparalleled archives of Phanerozoic Eon glacial oscillations (e.g., Miller et al., 2004, 2012; Isaacson et al., 2008; Rygel et al., 2008; Loi et al., 2010; Hoffman, 2011). However, tectonics, dynamic topography, and sedimentation can distort the preservation of depositional sequences, the apparent water-depth change of facies juxtapositions, and the number of erosional surfaces, which can lead to discrepant reconstructions of the number, magnitude, and relative age of inferred glacial cycles from continental shelf stratigraphy (Jervey, 1988; Einsele, 1993; Catuneanu et al., 2009). Many reconstructions of pre-Cenozoic glaciations that account for solid earth deformation and sedimentation presume that glacial oscillations impart uniform changes in accommodation and, therefore, identical stratigraphic architectures (Ghienne et al., 2014). Yet, a rich geophysical literature indicates that glacial isostatic adjustment (GIA), defined as the gravitational changes and viscoelastic deformation of Earth's crust due to the loading and unloading of ice sheets and oceans, generates spatially variable sea-level change and accommodation (Farrell and Clark, 1976; Milne and Mitrovica, 2008; Tamisiea and Mitrovica, 2011; Creveling et al., 2018).

Thus, shelf stratigraphy during glacial epochs represents a combination of the relative rates of (1) local sea-level (LSL) change imparted by GIA, (2) sedimentation, and (3) solid earth deformation. In this study, we held factors 2 and 3 constant to focus on the extent to which factor 1 controlled the local stratigraphic expression of the number, scale, and age of “glacioeustatic” oscillations and, hence, obscured the along-strike correlation of glacial-interglacial depositional sequences.

To explore the contribution of GIA to stratigraphic architectures, we developed a numerical model that creates computer-generated stratigraphic columns by combining LSL histories with identical model sedimentation histories under the assumption of no tectonic deformation. LSL histories refer to time series of local sea level, relative to the present, that capture spatial variability in GIA.

We computed the marine oxygen isotope stage (MIS) 11 to present-day (400 k.y.) LSL histories for nine bathymetric depth sites—from 1 to 250 m depth—along each of four shelf-perpendicular transects spanning the Pacific coast of North America, a region variably impacted by GIA during Quaternary glacial-interglacial cycles (Figs. 1A–1C; Fig. S1 in the Supplemental Material¹). These LSL predictions were generated from a gravitationally self-consistent sea-level theory (Kendall et al., 2005) that adopted the global ice-load history of Raymo and Mitrovica (2012) and the best-fitting regional one-dimensional (1-D) viscosity structure inferred by Creveling et al. (2017) for MIS 5 sea-level indicators. The model included an elastic lithosphere thickness of 95 km and uniform upper- and lower-mantle viscosities of 10²⁰ Pa·s and 5 × 10²¹ Pa·s, respectively. The LSL histories represent model estimates of local accommodation change, i.e., the change in sea surface versus radial surface of the crust due to GIA (Mitrovica and Milne, 2003) at 1 k.y. time steps, which we linearly interpolated to annual resolution. Sediment loading (Dalca et al., 2013) and compaction (Ferrier et al., 2017) were not incorporated.

We paired each of the 36 LSL histories (i.e., nine sites along four transects; Fig. S1) with the same 50 sedimentation histories. Each history drew annual depositional and erosional events (n = 400,000) from a truncated double-Pareto probability distribution function (PDF; Fig. 1D) for Cenozoic siliciclastic shelf sedimentation, with a mean of 30 ± 4 cm/k.y. (1σ; Trampush and Hajek, 2017). Model accumulation was comparable to late Pleistocene cores from Integrated Ocean Drilling Program (IODP) shelf sites within the area (Lyle et al., 2000; Akiba et al., 2009).

The numerical model generated synthetic stratigraphic columns on top of an initial bathymetry (set to modern) by tabulating 400 k.y. sediment accumulation at each site assuming no lithification (Fig. 2). For every model year, local accommodation, or the modeled water depth, dictated sediment accumulation (Fig. 2A; Jervey, 1988). Sediment accumulated, or eroded, when positive accommodation exceeded the thickness of an annual depositional or erosional sedimentation event (Fig. 2B). When a depositional event filled accommodation, any remainder of the event was discarded (Fig. 2C). Forced regression induced erosion to the elevation of the LSL lowstand (Fig. 2D). If subsequent LSL rise generated positive accommodation again, then the process was repeated (Fig. 2E). To facilitate the correlation of synthetic stratigraphies within individual transects, and laterally between transects, the model identified subaerial unconformities (SUs; where forced regression resulted in subaerial exposure), correlative conformities (CCs; where forced regression did not result in subaerial exposure), and maximum flooding surfaces (MFSs; where maximum accommodation occurred between two SUs or their CCs; Fig. 2F; Johnson and Murphy, 1984; Van Wagoner et al., 1987, 1988; Galloway, 1989; Hunt and Tucker, 1992; Posamentier and Allen, 1999). Only the preserved depositional events and surfaces were stacked into synthetic stratigraphic columns (Fig. 2G).

To quantitatively compare records, for each column, the model computed (1) a count of genetic stratigraphic sequences (GSSs), consisting of transgressive-regressive systems tracts (Fig. 2G; Galloway, 1989); (2) the difference in water depth between the oldest and youngest depositional event within each systems tract (Fig. 2G); and (3) the year of transgression after a SU within a given transect. By excluding physical processes such as loading, compaction, and lithification, and by asserting identical sedimentation histories and an initial topography, we did not seek to predict the real lithostratigraphy at these sites, but rather isolate the systematic effects of GIA on a local stratigraphic architecture. Incorporating local estimates for these features would enhance the geographic complexity imparted by GIA.

During glaciations, the imposed LSL history differed both within individual transects (Fig. S1) and between transects down the coastline (Fig. 1C). The loading of a growing ice sheet caused a depression of the underlying crust, displacing mantle material laterally to create a peripheral bulge (Fig. 1B; Farrell and Clark, 1976). Across a deglaciation phase, the subsequent melting of the ice sheet led to (1) postglacial rebound of the crustal depression and subsidence of the peripheral bulge, resulting in LSL rise, and (2) gravitational effects that caused water to migrate away from the ice sheet, culminating in LSL fall. Because factor 1 tended to dominate over factor 2, and the magnitude of both effects dropped with distance from the ice load, the net result was that deglacial LSL rise decreased southward across the subsiding outer flank of the peripheral bulge (i.e., transects 2–4 [T_2–4]). Additionally, the trend flattened out by the latitude of T₃, so the LSL predictions for T_3–4 showed reduced variability (Fig. 1C).

T₁ was located close to the edge of the ice load at maximum extent, and a large drop in LSL associated with gravitational effects combined with a more moderate crustal displacement across the deglaciation phase to culminate in a muted LSL change relative to T_2–4. The proximity of T₁ to the ice load introduced large shore-perpendicular gradients in the predicted sea-level signal (Fig. S1A). In contrast, peripheral bulge and gravitational effects at transects further south tended to have significantly smaller shore-perpendicular gradients, although ocean loading effects (Nakada and Lambeck, 1989; Mitrovica and Milne, 2003) may have contributed to these gradients (Figs. S1B–S1D).

Figure 3 illustrates a selection of the sequence-stratigraphic and lithostratigraphic architectures that arose from the pairing of MIS 11 to present-day LSL histories with one random, but representative, sedimentation history. From each depth-site library (n = 50), a comparison of synthetic stratigraphic columns within individual transects revealed that, as bathymetric accommodation increased downdip toward the basin, so too did overall stratigraphic thickness, the number of preserved GSSs (Figs. 3A and 3B), and the water-depth change within a regressive system tract (Figs. 3C and 3D).

A comparison of intermediate-depth (100–200 m) synthetic stratigraphic columns between transects (along depositional strike) exposed significant differences in these factors with a consistent geographic pattern (Figs. 3A and 3D). For example, T₁ consistently preserved the longest columns with the highest count of GSSs (e.g., four sequences at 130 m depth), which pinched out at T₂ to two sequences and expanded again toward T_3–4 to approximately three sequences (Fig. 3A). For all sedimentation scenarios, the difference in mean count of GSSs between transects was highest and statistically significant (p < 0.05) for the 130–200 m depth sites (Fig. 3B; Table S1). Similarly, between transects, the apparent change in water depth within a GSS varied widely (Fig. 3C); for example, T₂ recorded the largest local change in water depth leading into the Last Glacial Maximum (LGM) lowstand for sites 200 m and deeper and, thus, the largest apparent magnitude of ice-volume growth (Figs. 3C and 3D; Table S2). However, these geographic patterns in the number and magnitude of glacial cycles did not hold true for the 1–50 m and 250 m depth sites, as these columns were more statistically similar to each other (Figs. 3B and 3C), due to the limited and unlimited bathymetric accommodation, respectively, relative to the magnitudes of LSL variations (Fig. 1C).

Figure 4 presents actual age distributions for the year of preserved transgressive deposition atop the local MIS 2 (LGM) lowstand SU or CC at selected depth sites. Model runs revealed up to a 5 k.y., or ¹/₂₄ of the most recent glacial-interglacial cycle, updip diachroneity as sea-level rise reached progressively shallower depth sites (Fig. 4, same color distribution across depth sites). Across-shelf diachroneity was more pronounced (Fig. 4, different color distributions for the same depth site). For example, across the 130 m depth sites, local transgression at T₄ resumed ~8 k.y., or ¹/₁₅ of the cycle, earlier than at T₁ (16.89 ka ± 4.41 k.y. [1σ] vs. 8.95 ka ± 7 k.y., respectively). Wherever diachroneity existed for MIS 10, 8, and 6 transgressive surfaces, the local expression of a given surface could appear up to ~150 k.y., or >1 cycle, earlier at T₁ (Fig. S2).

If glacioeustasy were a global signal, and tectonics and sedimentation were everywhere identical, then stratigraphic columns between transects at the same bathymetric depths would be indistinguishable and retain the same information about sea-level and ice-volume change. Instead, intermediate-depth sites on the continental shelf displayed the largest departures from one another. This is a direct reflection of the contribution of GIA to the preservation and character of glacial-interglacial stratigraphic cycles produced by this forward model.

The library of stratigraphic records at the most ice-proximal intermediate site, T₁ at 130 m depth, consistently preserved the highest number of glacial cycles, the earliest ages for deglaciations (save for LGM), and, as expected from the LSL curve (Fig. 1C), the smallest magnitude of LSL change. In contrast, T₂ at 130 m depth, located on the outer flank of the peripheral bulge, recorded the least number of glacial cycles, and the magnitude of LSL change leading into the LGM was obscured by erosion (Figs. 3C and 3D). The intermediate-depth sites better revealed the retrogradation and progradation of lithofacies in response to GIA-induced LSL change at individual transects. However, even when sedimentation and tectonics were held identical across the coastline, GIA caused locations between transects to preserve a different number of glacial-interglacial cycles, apparent magnitude of LSL lowstand, and age of preserved transgressive deposition after a SU/CC in the stratigraphic record.

These model experiments reveal that GIA can complicate the along-strike correlation of continental shelf stratigraphic records across 10²–10³ km distances (T₁–T₂), and farther (T₃–T₄; Figs. 3A and 3D). Depending on the preserved field margin relative to the contemporaneous ice margin, a stratigrapher may conclude conflicting information about the same glacial epoch.

GIA creates geographically variable LSL change through glacial-interglacial cycles. In this study, we demonstrated how GIA-induced changes in LSL and accommodation alter the preservation of glacial-interglacial GSSs in continental shelf stratigraphic records. When tectonics, sedimentation, and bathymetric depths are identical, synthetic stratigraphic columns between transects display different information about the number, magnitude, and timing of ice-volume change. As a result, accurate correlation of continental shelf synthetic stratigraphic columns over several glacial oscillations is most difficult between sites proximal to the former ice sheet (T₁) and those on the outer flank of the peripheral bulge, though closest to the bulge crest (T₂). The accurate correlation of columns trending southward across the more distant sections of the outer flank of the peripheral bulge (T₃ and T₄) is easier. The range in bathymetric depths over which this pattern is observed will shift depending on initial bathymetry, sediment accumulation rates, and LSL histories. Improved understanding of the signal of GIA in continental stratigraphic records can refine reconstructions of the number, magnitude, and age of glaciations inferred from this physical archive.

This research was made possible by U.S. National Science Foundation award 2046244, a Geological Society of America (GSA) graduate student research grant, and the Oregon State University George and Danielle Sharp Fellowship. We thank J.X. Mitrovica for LSL histories and thoughtful discussion, and Steven Holland and two anonymous reviewers whose comments improved this paper.