The effect of lateral variations in Earth structure on Last Interglacial sea level

SUMMARY

It is generally agreed that the Last Interglacial (LIG; ∼130–115 ka) was a time when global average temperatures and global mean sea level were higher than they are today. However, the exact timing, magnitude and spatial pattern of ice melt is much debated. One difficulty in extracting past global mean sea level from local observations is that their elevations need to be corrected for glacial isostatic adjustment (GIA), which requires knowledge of Earth’s internal viscoelastic structure. While this structure is generally assumed to be radially symmetric, evidence from seismology, geodynamics and mineral physics indicates that large lateral variations in viscosity exist within the mantle. In this study, we construct a new model of Earth’s internal structure by converting shear wave speed into viscosity using parametrizations from mineral physics experiments and geodynamic constraints on Earth’s thermal structure. We use this 3-D Earth structure, which includes both variations in lithospheric thickness and lateral variations in viscosity, to calculate the first 3-D GIA prediction for LIG sea level. We find that the difference between predictions with and without lateral Earth structure can be metres to 10s of metres in the near field of former ice sheets, and up to a few metres in their far field. We demonstrate how forebulge dynamics and continental levering are affected by laterally varying Earth structure, with a particular focus on those sites with prominent LIG sea level records. Results from four 3-D GIA calculations show that accounting for lateral structure can act to increase local sea level by up to ∼1.5 m at the Seychelles and minimally decrease it in Western Australia. We acknowledge that this result is only based on a few simulations, but if robust, this shift brings estimates of global mean sea level from these two sites into closer agreement with each other. We further demonstrate that simulations with a suitable radial viscosity profile can be used to locally approximate the 3-D GIA result, but that these radial profiles cannot be found by simply averaging viscosity below the sea level indicator site.

Composition and structure of the mantle, Mantle processes, Sea level change, Dynamics of lithosphere and mantle, Rheology: mantle

1 INTRODUCTION

The Last Interglacial (LIG; ∼130–115 ka) is a time in Earth’s history during which global average temperatures were 1–2 °C warmer than pre-industrial values (Dutton et al. 2015a). As such, it has been used as a testing ground to study how ice sheets and sea level respond to past and possibly future warming (DeConto & Pollard 2016; Fischer et al. 2018). Reconstructions of global mean sea level (GMSL) during the LIG are based on sea level indicators, such as fossil corals, that constrain the local elevation of sea level at their time of formation (Rovere et al. 2016). Once locally reconstructed, this elevation has to be corrected for processes that cause a deviation between local sea level and the GMSL. One of these processes is glacial isostatic adjustment (GIA), which is the response of Earth’s solid surface, gravity field and rotation axis to changes in ice and ocean mass. GIA is an important contributor to interglacial sea level change, even far away from major ice sheets (e.g. Mitrovica & Milne 2002; Lambeck et al. 2012). In addition to GIA, other processes such as earthquakes, crustal deformation, sediment loading and dynamic topography can further deform Earth’s surface and cause local sea level change (Briggs et al. 2014; Austermann et al. 2017; Pico 2019; Stephenson et al. 2019).

Estimates of peak GMSL during the LIG, whether based on records from individual sites (O’Leary et al. 2013; Dutton et al. 2015b) or by combining data from multiple locations into a statistical framework (Kopp et al. 2009), are generally between 6 and 9 m. However, some recent work suggests that this range may be overestimating GMSL during the LIG (Clark et al. 2020; Dyer et al. 2021). In general, significant debate continues about both the magnitude of excess melt (relative to present-day) and its timing: data from Western Australia indicate that GMSL exceeded present-day values by a few metres early in the LIG, followed by a GMSL rise up to 9 m towards the end of the LIG (O’Leary et al. 2013). This reconstruction is in disagreement with estimates obtained from the Seychelles, which indicate that high GMSL was attained early in the LIG and continued to slowly increase, with possible intermittent sea level drops (Dutton et al. 2015b; Vyverberg et al. 2018). While constraints from late LIG sea level are absent in the Seychelles, they are present in Xcaret, Mexico—where sea level has been argued to undergo a step increase around 121 ka (Blanchon et al. 2009)—and Mallorca, where speleothem records indicate constant or slightly falling GMSL throughout the LIG (Polyak et al. 2018). A global compilation of data indicates an oscillation in sea level with a highstand both early and late (Kopp et al. 2009), however there's a lack of evidence for this evolution in proximal ice records and ice sheet dynamics (Barlow et al. 2018).

Ongoing disagreement regarding the magnitude, timing and spatial distribution of LIG melt raises the possibility that complexities associated with the GIA correction may be responsible for some of these differences. As noted above, each local sea level estimate needs to be corrected for GIA to infer GMSL. The GIA correction requires both an ice history and a viscoelastic structure for Earth’s interior as input, both of which are underconstrained. Uncertainties associated with the ice history can change the GIA correction by several metres during the LIG (Lambeck et al. 2012; Dendy et al. 2017; Rohling et al. 2017). In regard to Earth’s viscoelastic structure, previous studies of LIG sea level have all assumed that Earth’s viscosity varies purely as a function of depth. However, based on evidence from seismic tomography, mineral physics and geodynamics, it is expected that significant lateral variations exist in both Earth’s viscosity and lithospheric structure (e.g. Dannberg et al. 2017; Priestley et al. 2018). Indeed, these lateral variations are important for understanding the impact of GIA on sea level during the last deglaciation (Austermann et al. 2013; Li et al. 2018; Kuchar et al. 2019) and affect the pattern of present-day deformation across Antarctica (Nield et al. 2018; Gomez et al. 2018).

In this study, we focus on investigating how lateral variations in Earth structure affect sea level during the LIG. We generate a new model of lateral Earth structure that is based on seismic tomography. In contrast to previous work, which adopts a pre-determined scaling from shear wave speed into viscosity for the upper mantle (Austermann et al. 2013; Gomez et al. 2018; Li et al. 2018), we invert laboratory-based parametreizations for material properties using a suite of independent constraints on mantle structure (Richards et al. 2020). We pair this earth model with an ice history to predict the effect of lateral viscosity variations on sea level at key sites, and provide physical insights into the GIA changes predicted both in the near and far field (i.e. close and distant to the ice sheet). Given the computational expense of such calculations, we are limited to performing a relatively small set of exploratory simulations, but these nevertheless provide a first estimate of the potential magnitude and geometry of the LIG GIA signal associated with realistic departures from radial mantle viscosity profiles. While this work is focused on the LIG, insights on the physical mechanisms hold for any interglacial period and are therefore also relevant to earlier interglacials (e.g. MIS 11) and the late Holocene. We also investigate how well the GIA signal obtained when including lateral variability in Earth structure can, at a given location, be accurately represented by a GIA model with a purely radial Earth structure. Finally, we compare our results to LIG sea level records at key sites to consider the extent to which lateral Earth structure and the associated GIA prediction might impact estimates of GMSL over the LIG.

2 METHODS AND DATA

2.1 GIA numerical model

To investigate how GIA causes spatially varying sea level over the LIG, we use a finite volume-based approach to solve for radial displacement of the solid Earth, its change in gravity field, and rotation axis in response to an evolving ice and ocean load (Latychev et al. 2005). The numerical approach incorporates lateral variations in Earth structure and calculates the resulting gravitationally self-consistent sea level change, while accurately accounting for shoreline migration (Mitrovica & Milne 2003). This computational model is well established, having been used in many previous studies (e.g. Austermann et al. 2013; Goldberg et al. 2016; Gomez et al. 2018; Kuchar et al. 2019). GIA calculations described here for radially symmetric Earth structure are performed using both the finite volume approach described above and the pseudo-spectral approach described in Kendall et al. (2005).

2.2 Ice reconstruction

The ice reconstruction we adopt is based on a combination of several published studies in order to obtain satisfactory fits to multiple independent data sets (Fig. 1). From the present-day back to 26 ka, we use the deglacial ice sheet reconstruction ICE-6G (Peltier et al. 2015). For the preceding glaciation, we use the reconstructions by Pico et al. (2017) and Creveling et al. (2017), which are more consistent with sea level observations from these time periods. To isolate the GIA contribution to sea level during the LIG, we assume present-day ice geometry from 128 to 117 ka. Any excess ice melt beyond the present-day level will, of course, produce an additional geographically varying sea level fingerprint (Hay et al. 2014). The timing and melt geometry of the penultimate deglaciation has been widely debated, with estimates including a smaller (Rohling et al. 2017) or larger (Shakun et al. 2015) overall ice volume compared to the last deglaciation; a deglaciation that commenced early (around 140 ka; Thomas et al. 2009) or late (closer to 135 ka Waelbroeck et al. 2002); and an ice distribution characterized by a slightly (Lambeck et al. 2006) or significantly (Colleoni et al. 2016) larger Fennoscandian ice sheet relative to LGM. The ice geometry and timing of melt across the penultimate deglaciation will significantly affect GIA during the LIG and this has been explored in detail elsewhere (Dendy et al. 2017; Rohling et al. 2017). Here, we adopt a representative ice history that has a total ice volume at the penultimate glacial maximum (PGM) that is similar to the last glacial maximum (LGM; consistent with oxygen isotope estimates Waelbroeck et al. 2002), a slightly slower deglaciation (in line with U–Th dated corals from Tahiti Thomas et al. 2009), and an ice distribution characterized by a larger Fennoscandian ice sheet, following Lambeck et al. (2006), and in turn smaller Laurentide ice sheet (Fig. 1). Our calculations start at 150 ka and continue until the present-day.

Figure 1.

Ice model reconstruction. (a) Distribution of northern hemisphere ice during the penultimate glacial maximum (PGM) used in our model. The Fennoscandian ice sheet is based on Lambeck et al. (2006). B) Global mean sea level with respect to present-day over the course of the model. Vertical coloured bars indicate timings of the last glacial maximum (LGM), Last Interglacial (LIG) and penultimate glacial maximum.

2.3 Elastic Earth structure

To model the instantaneous elastic deformation of the solid Earth in response to evolving ice loads, we require estimates of the bulk and shear moduli of the mantle. We adopt the radially symmetric model STW105 (Kustowski et al. 2008), which uses more data and an improved inversion scheme and crustal correction compared to earlier studies such as PREM (Dziewonski & Anderson 1981). We emphasize that the choice of the elastic structure has only a minimal effect on our results (i.e. on the order of centimetres during the LIG). While elastic moduli are known to vary laterally, these perturbations are considerably smaller than those in viscosity and do not play a first-order role in the GIA response. We therefore maintain a radially symmetric elastic structure and only investigate the effect of lateral viscosity variations in this study.

2.4 Viscous Earth structure including lithospheric thickness variations

The mantle convects with a time-dependent planform that evolves on million-year timescales and introduces significant lateral temperature heterogeneity (Turcotte & Schubert 2002). Laboratory experiments on the deformation of mantle rocks show that their viscosity is strongly dependent on temperature, and it has therefore long been known that significant lateral variations in viscosity occur within the mantle (Cathles 1975; Ranalli 1995). The velocity of seismic waves is also sensitive to mantle temperature and rheology, and seismic tomography is therefore our most direct tool for imaging mantle structure (Bullen 1975).

Seismic velocities are traditionally converted into temperature and then viscosity using a combination of physical and phenomenological laws (e.g. thermal expansion, an Arrhenius relationship to describe the temperature dependence of viscosity) and material properties that have been derived from laboratory experiments. Here, we construct a mantle viscosity structure that relies on these same conversion relationships, including an up-to-date treatment of anelasticity at seismic frequencies that is detailed below. As with previous studies, uncertainties in the appropriate material properties, the rheological deformation mechanism responsible for accommodating GIA motions, and variations in measured seismic velocity structure between different tomography models all introduce uncertainty into the resulting viscosity structure. In contrast to other studies, however, we leverage additional information on the thermal and rheological state of the upper mantle to optimize the selection of appropriate material properties. This approach substantially reduces uncertainties in inferred mantle viscosity structure, which is demonstrated and discussed in Section 4.1.

2.4.1 Viscosity above the transition zone and lithospheric thickness

When a polycrystalline viscoelastic material - such as the mantle - is cold, deformation associated with the passage of acoustic energy is elastic, yielding a linear dependence of shear wave velocity (V_S) on temperature referred to as the anharmonic velocity. As temperature increases, however, anelastic deformation (a special case of fully recoverable viscoelastic deformation) also begins to occur due to the presence of point defects, dislocations and grain boundaries. This additional process results in a non-linear relationship between V_S and temperature and is particularly important to account for when inferring viscosity in high-temperature regions (Karato 1993; Wu et al. 2012). Anelastic behaviour has been extensively studied in laboratory experiments on silicates and organic analogues of mantle rocks, revealing that the strength of the anelastic regime varies with both the frequency of seismic waves and as a function of material properties, such as melting temperature and grain size (Sundberg & Cooper 2010; McCarthy et al. 2011; Faul & Jackson 2015). Several studies have attempted to parametreise these complex dependencies and have been regularly updated as forced oscillation and creep experiments in the laboratory have been pushed towards increasingly realistic frequencies, pressures, temperatures, grain sizes and strain rates (Jackson & Faul 2010; Takei 2017).

In this study, we map V_S into temperature and viscosity in the upper 400 km of the mantle using the parametreization of Yamauchi & Takei (2016), which has been developed through forced oscillation experiments on borneol. The parametreization accounts for the effects of anelasticity in pre-melt conditions, when temperature (T) exceeds ∼90 per cent of the melting temperature (T_m; both defined in Kelvin). These conditions most often occur in regions of the asthenosphere that underlie thin lithosphere, such as beneath West Antarctica, which is a site that experiences significant ice mass changes over the glacial cycle. Specific details on the representation of anelasticity are provided in Appendix A. To summarize, seismic velocity and attenuation are self-consistently tied to temperature and steady-state diffusion creep viscosity via a system of coupled equations that depend on seven material properties (including the activation energy, E_a, which controls the dependence of viscosity on temperature through an Arrhenius relationship). Here, we assume that temperature is the dominant cause of seismic velocity variations and that grain size and composition play only a second-order role. The equations that are being used in this study are given by eqs (3)–(17) of Richards et al. (2020).

The standard approach is to adopt material properties and their associated uncertainties that are appropriate for upper mantle rocks (typically olivine) and have been obtained from laboratory experiments (e.g. Kaufmann et al. 2005; van der Wal et al. 2010; Li et al. 2018). Here, rather than fixing these properties using an assumed mineralogy, we take advantage of an inverse calibration scheme outlined in Richards et al. (2020) that considerably reduces uncertainty in inferred mantle structure (see Section 4.1). The philosophy behind the approach is that certain physical properties of the Earth are ‘known’, including the typical thermal structure of oceanic lithosphere, the average adiabatic temperature profile within the convecting mantle, the attenuation structure of the upper mantle beneath old oceanic lithosphere, and the mean diffusion creep viscosity of the upper mantle from studies of GIA. Any model of upper mantle temperature and viscosity structure inferred from shear wave velocities should be compatible with these constraints, and we therefore restrict ourselves to a subset of material properties that also satisfy these physical characteristics.

To generate the constraints, we first stack shear-wave velocities from the tomography model (described below) as a function of depth and oceanic age. Temperature contours from the Richards et al. (2018) plate cooling model are subsequently overlain and V_S–T tie points at depths of 75, 100 and 125 km are extracted. A second set of deeper tie points is generated by assuming that the average value of V_S as a function of depth over the 225–400 km range should yield a temperature that is consistent with the 1333 °C adiabat. A third constraint is obtained by overlaying the observed attenuation structure at depths of 150–400 km beneath >100 Ma seafloor from the QRFSI12 model of Dalton et al. (2009) on the equivalent V_S stack, in order to generate a set of V_S–Q⁻¹ tie points as a function of depth. Finally, we require that the mean diffusion creep viscosity from 225 to 400 km depth be equivalent to the average upper mantle value of 3 × 10²⁰ Pa s that has been obtained from previous studies of GIA (Lau et al. 2016). We calculate a range of predicted temperature, attenuation and viscosity maps by varying the seven material properties and comparing the results to the calibration targets described above. Misfit is optimized by iteratively updating the material properties. To reduce the likelihood of locating local minima in the inversion, we use a two-step minimization algorithm consisting of an initial, relatively coarse parameter sweep followed by Powell’s conjugate gradient algorithm. Further details on this calibration scheme can be found in Richards et al. (2020).

Here, we use the SL2013sv tomography model, which has global coverage in the upper mantle, and patch in the SL2013NA regional update in North America that takes advantage of the dense station coverage afforded by the USArray seismic network (Schaeffer & Lebedev 2013, 2014). We have chosen this model for two reasons. First, it has high horizontal resolution (∼280 km horizontal node spacing) and is constructed from both body waves and a large quantity of surface wave data (including higher modes), which are particularly sensitive to velocity structure in the upper ∼350 km of the mantle. Secondly, density and temperature fields derived from this model have been shown to correlate well with independent geophysical and geological observations including gravity anomalies, residual topography, continental geotherms and mineral deposits (Steinberger 2016; Hoggard et al. 2017, 2020). The anelastic calibration scheme yields optimal material properties of 74.7 GPa for the reference shear modulus (with 1σ uncertainties of approximately 3 per cent). Its dependence on temperature is –16.1 MPa °C⁻¹ (∼12 per cent) and on pressure is 2.56 (dimensionless, ∼7 per cent), respectively. The reference diffusion creep viscosity is 2.51 × 10²¹ Pa s (∼8 per cent), and its dependence on temperature and pressure are controlled by an activation energy of 304 kJ mol⁻¹ (±250 kJ mol⁻¹) and activation volume of 3.0 cm³ mol⁻¹ (±6.0 cm³ mol⁻¹). All uncertainties are 1σ. We note that a negative activation volume would imply that the sensitivity of viscosity to temperature may decrease with depth, which has previously been suggested for mantle mineral assemblages when self-diffusion of certain ions (such as Si and O) becomes rate-limiting (Fei et al. 2018; Jain et al. 2019). The solidus gradient is 0.946 °C km⁻¹ (∼25 per cent). The resulting scaling relationship between shear wave velocity and lateral viscosity perturbations is shown in Fig. 2. Throughout the manuscript we define lateral viscosity perturbations as |$\log_{10}\frac{\eta _1}{\eta _0}$|⁠. The effect of accounting for anelastic effects in this manner is to increase viscosities by between 1 and 1.5 orders of magnitude in the slowest V_S, lowest-viscosity locations, in comparison to a purely anharmonic conversion (Fig. 2e). Viscosities are unaffected in faster velocity regions where temperatures are colder (i.e. T < 0.9T_m).

Figure 2.

Conversion of shear wave speed to viscosity perturbations. (a) Calibrated relationship between shear wave velocity (V_S) and lateral viscosity perturbations as a function of depth for our calibrated upper-mantle model of anelasticity at seismic frequencies (Schaeffer & Lebedev 2013, 2014; Yamauchi & Takei 2016; Richards et al. 2020). Coloured lines transition to dashed grey at the lithosphere-asthenosphere boundary (i.e. 1175 °C); black dotted lines are isothermal contours from 1200 to 1600 °C at 50 °C increments. (b) Same if anelastic effects are excluded. (c and d) Same values as in (a) and (b), respectively, but plotted as V_S versus depth and coloured by viscosity perturbations. (e) Difference between (c) and (d), illustrating that including anelastic effects primarily acts to increase the inferred viscosity at slow shear wave velocities. Throughout this study, we include anelastic effects and use the conversion shown in panels (a) and (c).

Lithospheric thickness is identified using the depth of the 1175 °C isothermal surface, which has an average global value of ∼100 km and varies from ∼300 km in cratonic regions to <25 km in sites of active rifting and at mid-oceanic spreading centres (Fig. 3a, Hoggard et al. 2020). This specific isotherm is chosen because it coincides with seismological evidence for the depth of the lithosphere–asthenosphere boundary in oceanic regions (Richards et al. 2018). At 100 km depth, a temperature of 1175 °C corresponds to a diffusion creep viscosity of 4.5 × 10²² Pa s in our anelasticity parametreization, which yields a characteristic Maxwell time of ∼20 kyr. We therefore consider this boundary an adequate representation of the transition between asthenospheric material that undergoes viscoelastic deformation during the simulation, and lithospheric material that only deforms elastically. It has been suggested that lithospheric thickness on GIA timescales appears thinner than inferred from seismic tomography due to the onset of viscous or transient deformation (Watts et al. 2013; Lau et al. 2020). We therefore perform an additional sensitivity test where lithospheric thickness is reduced everywhere by 20 per cent (i.e. 80 km global average).

Figure 3.

3-D viscosity structure of the Earth. (a) Global lithospheric thickness variations, yielding a globally averaged lithospheric thickness of 100 km (colourbar at the bottom left). (b–f) Lateral viscosity perturbations at depths of 200, 300, 400, 600 and 1200 km, respectively (colourbar at the bottom right). Thick black contour in panel (b) delineates lithospheric portions.

2.4.2 Viscosity within and beneath the transition zone

Deeper than 400 km, the sensitivity of surface waves to velocity structure drops significantly, the dominant mantle mineralogy switches from olivine to wadsleyite, and the number of independent constraints on mantle properties is considerably more limited. For these reasons, we use a different approach to estimate viscosity perturbations below 400 km depths. For the tomography model, we take the whole mantle SEMUCB-WM1 model, which is constructed using a combination of body and surface wave data and uses a sophisticated hybrid approach to numerically simulate wavefield propagation and invert waveforms for shear wave-velocity structure (French & Romanowicz 2014). We linearly blend the upper and lower shear wave model over the 300–400 km depth range.

To generate lower mantle viscosity perturbations, we first convert V_S into temperature using the Perple_X Gibbs free-energy minimization software of Connolly (2005) and the thermodynamic database of Stixrude & Lithgow-Bertelloni (2011). Assuming a pyrolitic composition, the mineralogical make-up is estimated as a function of pressure and temperature, yielding predictions of the elastic moduli and density that can be converted into anharmonic velocity. At a given depth and for an individual mineral assemblage, these material properties exhibit an approximately linear dependence on temperature, resulting in a unique conversion from anharmonic velocity into temperature. In the vicinity of phase transitions, however, a discontinuity occurs that can render this conversion non-unique. To avoid this issue, at each depth, we first linearize the temperature-dependence of the density and elastic moduli over a temperature window that extends ±500 °C around a quasi-steady state geotherm obtained from mantle convection simulations (Supplementary Material, Schuberth & Bunge 2009). The anharmonic velocity as a function of pressure and temperature is subsequently corrected for anelastic effects using the Q5 radial attenuation profile, eqs (1) and (2) of Cammarano et al. (2003), and a mantle solidus from Andrault et al. (2011). To ensure that we obtain a self-consistent mantle geotherm from the tomography model, we extract the V_S variation over the ±500 °C temperature window and pin the absolute V_S at the geotherm temperature to the average V_S of the tomography model at this depth. Thus lateral V_S variations at any depth are mapped into temperature variations around the average geotherm. Finally, these temperature variations are converted into viscosity perturbations using a radial activation enthalpy profile constructed from our upper mantle result extended down to 660 km, followed by the lower mantle profile of Steinberger & Calderwood (2006).

2.4.3 Resulting viscosity structure

Our resulting lithospheric thickness and viscosity perturbations at specific depths are shown in Fig. 3 and are provided as supplementary material. As expected, thick lithosphere aligns with cratonic regions and it thins towards mid-ocean ridges (Fig. 3a). Oceans are generally underlain by a less viscous asthenosphere and higher viscosity slabs appear beneath subduction zones at greater depth. The general patterns and order of magnitude viscosity variations are comparable to previously derived viscosity structures (Li et al. 2018). For our GIA calculations, we superimpose the lateral viscosity perturbations shown in Fig. 3 on different 1-D viscosity profiles. Our default simulations use a 1-D viscosity profile referred to as p55, which consists of an upper mantle viscosity of 5 × 10²⁰ Pa s and a lower mantle viscosity of 5 × 10²¹ Pa s (Raymo et al. 2011). Uncertainties in inferred Earth structure and their impacts on GIA predictions are addressed in Section 4.

3 RESULTS

We performed GIA calculations using both radially symmetric and laterally varying viscosity structure, which we will refer to as 1-D and 3-D simulations, respectively. Results and predictions from these simulations will be referred to as 1-D and 3-D results and predictions. The 1-D viscosity model is identical to the average of the 3-D viscosity model at each depth (averaging is done in log space) except in the analysis described in Section 4.3. We first describe and explain the general patterns that are common to both 1-D and 3-D simulations. We next detail differences between 1-D and 3-D results in the near and far field, before comparing results at specific sites with important sea level indicators from the LIG.

3.1 General patterns of GIA over the LIG

Near field. GIA models predict that sea level changes over the course of an interglacial will vary significantly across the globe (Raymo & Mitrovica 2012; Dutton & Lambeck 2012; Dendy et al. 2017). At the beginning of the LIG (Fig. 4a), relative sea level (i.e. sea level relative to present-day, Mitrovica & Milne 2003) is high in formerly glaciated regions because these areas are experiencing ongoing postglacial rebound in response to the just-completed deglaciation, while the peripheral bulges surrounding them are subsiding from an elevated position, leading to low values of relative sea level. That is, the solid Earth beneath the major ice sheets is in greater isostatic disequilibrium at the beginning of the LIG compared to today. For the Laurentide ice sheet, this pattern is reversed at the end of the LIG (Fig. 4b) at which time more rebound (and peripheral bulge subsidence) has occurred in comparison to today. This is the case because the LIG was longer than the present interglacial and because our ice history adopts a Laurentide ice sheet that was smaller during the PGM than the LGM. In contrast, the results for the Fennoscandian region suggest more isostatic disequilibrium at the end of the LIG compared to today (Fig. 4b)—the formerly glaciated area remains below present levels while the peripheral bulge sits above present levels. In this case, the greater size of the Fennoscandian ice sheet at PGM relative to the LGM more than compensates for the fact that the duration of the LIG was longer than the duration of the present interglacial.

Figure 4.

LIG sea level predictions assuming 1-D and 3-D Earth structure. (a and b) Prediction of relative sea level at the beginning (128 ka) and end (117 ka) of the Last Interglacial, respectively, assuming radially symmetric Earth structure. (c and d) Same as (a) and (b) but including lateral variations in lithospheric thickness and mantle viscosity. The a–a’ transect indicated in panel (a) is used in Fig. 6.

Far field. First-order sea level patterns in the far field include continental levering, ocean syphoning and rotational effects (Mitrovica & Milne 2002). Continental levering—a tilting of the crust at continental shorelines—is driven by the loading and unloading of the oceans by the water produced from ice melting, while the adjacent continents experience no such load changes. During interglacials, the process leads to a progressive rise in sea level on the oceanward side of a continental shoreline and a fall on the continent side of the shoreline. The length scale of this effect is related to the thickness of the lithosphere, with thicker lithosphere leading to a broader levering signal that extends further away from the shoreline. The levering process is significantly more advanced at present-day relative to the start of the LIG, leading to the strong gradient in sea level on a transect across most shorelines (Fig. 4a). In contrast, at the end of the LIG, the levering process is somewhat more advanced than at present-day, and the gradient along the transect is therefore of opposite sign and smaller in magnitude (Fig. 4b). In the middle of ocean basins, ocean-syphoning—the migration of water away from such regions and largely toward subsiding peripheral bulges—dominates sea level physics during interglacials and drives a sea level fall. At the beginning of the LIG, this process is less advanced than at present-day and sea level is higher than at present (Fig. 4a), while at the end of the (longer duration) LIG, the opposite is true and sea level is lower (Fig. 4b).

Models that include lateral variability in Earth structure exhibit the same first-order patterns described above (Figs 4c and d). However, there are notable shorter wavelength differences that arise due to lateral variations in both lithospheric structure (Figs 5a and b), and mantle viscosity (Figs 5c and d). The magnitude of the difference between 1-D and 3-D predictions is largest in the near field, where it is on the order of several to tens of metres. The difference is smaller, but still up to few metres, in the far field of ice sheets (Figs 5e and f). In the next section, we analyse in more detail the differing GIA responses.

Figure 5.

Effect of lateral variations in Earth structure on LIG sea level. (a and b) Difference in predicted sea level at the beginning (128 ka) and end (117 ka) of the Last Interglacial (LIG), respectively, between a simulation with variable lithospheric thickness above a radial viscosity profile and one with constant lithospheric thickness. (c and d) Difference in predicted sea level at the beginning (128 ka) and end (117 ka) of the LIG, respectively, between a simulation with variable sublithospheric viscosity variations and one with radially symmetric viscosity (both simulations have a constant lithospheric thickness). (e and f) Difference in predicted sea level at the beginning (128 ka) and end (117 ka) of the LIG, respectively, between a simulation including lateral variations in both lithospheric thickness and viscosity, and a purely radial model.

3.2 Near-field effects of lateral variations in Earth structure

Lithospheric thickness variations alone. Both regions that were covered by massive ice sheets during the PGM, Fennoscandia and Canada, are largely cratonic regions with a thick continental lithosphere (Fig. 3a). This similarity leads to similar patterns in North America and northern Europe when comparing the 1-D and 3-D results (Fig. 5a). The thicker lithosphere underneath the former Laurentide ice sheet leads to less subsidence and a more distal peripheral bulge during the glacial maxima. At the beginning of the LIG, when ice sheets were recently melted, this leads to higher topography (or lower sea level) in the centre of the former ice sheet with a broader peripheral bulge (Figs 5a and 6a). Towards the end of the LIG, as the peripheral bulge continues to subside, the difference between the 1-D and 3-D results is small (Figs 5b and 6b). Sea level in Fennoscandia exhibits a similar pattern at the beginning of the LIG except in western Europe (UK, Germany and France), which can be explained by the much thinner lithosphere in this region (Fig. 3a). The predicted response evolves towards the end of the LIG as the centre of rebound shifts slightly northwards (Fig. 5b).

Figure 6.

Near-field GIA effects due to lateral variations in Earth structure. (a and b) Relative sea level at the beginning (128 ka) and end (117 ka) of the LIG across the peripheral bulge of the Laurentide ice sheet (a–a’ transect of Fig. 4a). Results are shown for simulations using a 1-D Earth structure (black line) with a constant lithospheric thickness of 100 km, upper mantle viscosity 5 × 10²⁰ Pa s, and lower mantle viscosity 5 × 10²¹ Pa s. The dotted and dashed black lines show the prediction using lateral variations in lithosphere alone and full 3-D Earth structure (lithosphere plus mantle viscosity), respectively. (c and d) Relative sea level using a 1-D and 3-D simulations along the same transect at different times during the deglaciation. (e) Cross section of Earth structure along the transect from south to north. Viscosity perturbations are relative to the average upper mantle viscosity of 5 × 10²⁰ Pa s.

Lateral variability in both lithospheric thickness and mantle viscosity. We next consider the effects of including lateral variations in mantle viscosity in addition to lithospheric thickness variations. In Fennoscandia, the patterns change only slightly, most notably in western Europe (UK, Germany and France). These areas have a thin lithosphere but are underlain by mantle at 300–400 km depth that is 1–2 orders of magnitude more viscous than the global average. This high-viscosity feature, which is a slab associated with the closure of the Tethys Ocean (specifically the Apennine–Calabrian–Maghrebides slab, Fichtner et al. 2013; van Hinsbergen et al. 2014), causes a widening of the peripheral bulge observed in the 3-D to 1-D difference, both at the beginning and end of the LIG (Figs 5c–f). Considering the peripheral bulge of the Laurentide ice sheet along the U.S. East coast, we find that lateral viscosity variations bring the location of the bulge closer to the former ice sheet (Figs 5c and d) while remaining similar in amplitude (Fig. 6a). A significant difference between the 3-D simulation and the 1-D result is visible at the southern end of the peripheral bulge (i.e. south of ∼37 °N), which is underlain by low viscosity mantle at ∼300 km depth (Fig. 6e). This weak region has previously been associated with active upwelling flow above the Farallon slab (Rowley et al. 2013). We speculate that the low viscosity in this region focuses deformation associated with the peripheral bulge and possibly also continental levering, leading to a relative sea level high and more northern peripheral bulge in comparison to runs that only account for lithospheric thickness variations (Fig. 6). This leads to a relative sea level high that persists throughout the deglaciation and LIG (Figs 4c–f and 6a–d). This response is also visible along the U.S. west coast, which is generally underlain by lower viscosity mantle. It is, however, more localized at the edge of the former ice sheet, with a sea level peak occurring around 47 °N.

3.3 Far-field effects of lateral variations in Earth structure

Lithospheric thickness variations alone. Many far-field sea level sites are located at continental margins that sit on the transition from thinner oceanic lithosphere to thicker continental lithosphere (Fig. 3b). This differs for ocean island sites, which are generally situated on thinner oceanic lithosphere and can be underlain by lower viscosity if their origin is plume related. Lithospheric thickness variations affect far-field sea level indicators more at the beginning of the LIG than the end (Figs 5a and b) and in two specific ways: First, thicker lithosphere leads to continental levering over a broader area rather than being focused in a narrow corridor along the coastlines. When the lithosphere is thicker on the landward side of the coastline, the continental levering signal becomes asymmetric. One location that exemplifies this situation is Western Australia (Figs 7a and c). The amount of continental levering when including lateral variations in lithospheric thickness in the calculation is close to the 1-D prediction on the oceanward side, where the lithosphere is only slightly thicker than 100 km. However, predictions on the landward side are lower in amplitude and wider due to the thick (∼200 km) lithospheric root associated with the Yilgarn and Pilbara cratons (Fig. 7c). Secondly, large oceanic islands generally experience more continental levering when lateral variations in lithospheric structure are considered (Fig. 7d). This effect occurs because the lithosphere is typically thinner in oceanic settings than the 100 km global average. For example, the Seychelles are part of a granitic plateau in the western Indian Ocean that was exposed during the LGM and have a spatial extent of 50–100 km (Fig. 7b). The lithospheric thickness here is around 80–90 km, and undergoes continental levering during exposure. Thinning the lithosphere further causes the levering to become more pronounced, while increasing it to 100 km or higher results in the loss of continental levering effects (Fig. 7d, Dendy et al. 2017).

Figure 7.

Far-field GIA effects due to lateral variations in Earth structure. (a and b) Bathymetry around Australia and the Seychelles, respectively. Contour lines show lithospheric thickness in km. Relative sea level along the transects shown as black lines (b–b’ and c–c’) are displayed in panels (c) and (d), respectively. (c and d) Relative sea level is extracted at the beginning of the Last Interglacial (128 ka) and shown for simulations using a 1-D Earth structure (black line) with a lithospheric thickness of 100 km, upper mantle viscosity 5 × 10²⁰ Pa s, and lower mantle viscosity 5 × 10²¹ Pa s. The dotted and dashed black lines show predictions using lateral variations in lithosphere alone and full viscoelastic structure, respectively. The blue dashed line shows predictions using full 3-D Earth structure (lithosphere plus mantle viscosity), but with a thinner lithosphere that has a global average of 80 km instead of 100 km.

Lateral variability in both lithospheric thickness and mantle viscosity. Understanding the far-field response to full 3-D variability in Earth structure is challenging because far-field sea level indicators are not only sensitive to local Earth structure, but also to structure beneath the melting ice sheets and their periphery and to deeper mantle structure along the path between these ice sheets and the far-field site (Crawford et al. 2018). In Western Australia, we observe that including lateral viscosity variations leads to a smoother continental levering signal (Fig. 7c). The 1-D and 3-D simulations exhibit long-wavelength, metre-amplitude differences throughout the ocean basins, including a more positive sea level at the Seychelles (Fig. 7d). The slightly less pronounced continental levering in the full 3-D simulation is due to the higher than average viscosity beneath the Seychelles plateau. Reducing the global average lithospheric thickness to 80 km instead of 100 km (while keeping lateral variations in viscosity the same) allows for more deformation related to continental levering and leads to an increased sea level signal on the Seychelles plateau (Fig. 7d).

3.4 Predictions of 3-D GIA at key LIG sites

Next we consider predictions at specific sites that have notable records of LIG sea level (Fig. 8).

Figure 8.

GIA time series at key sites. (a–f) Local sea level at the Seychelles (La Digue), Western Australia (Cape Range), Yucatan Peninsula (Xcaret), Bahamas (Great Inagua and Eleuthera, Whale Point), Mallorca (Coves del Pirata) and Bermuda (Grape Bay), respectively. Thick grey line is the eustatic (global mean sea level) value assumed in the GIA model. Black lines show predictions of local sea level using a 1-D (solid line) and 3-D (dashed line, higher activation energy of 560 kJ mol⁻¹ in the asthenosphere—dotted line) GIA model on top of the p55 average radial viscosity profile. The yellow line in panel (d) shows the 3-D and 1-D GIA predictions for Eleuthera (Bahamas) instead of Great Inagua (Bahamas). Blue lines are the same as black lines but use an average lithospheric thickness of 80 km. Red lines are the same as black lines but use the background 1-D viscosity VM5 (Peltier et al. 2015). G) Locations from (a) to (f) shown on a map.

Near-field locations. Bermuda and Mallorca are located on the peripheral bulge of the former Laurentide and Fennoscandian ice sheets, respectively. This forebulge subsides over the course of the LIG and therefore leads to sea level rise if GMSL is assumed to be constant (solid black lines, Figs 8e and f). Accounting for lateral variations in viscosity at Mallorca leads to a larger sea level rise over the LIG (dashed black line, Fig. 8e). Bermuda, on the other hand, is located in a region that is not strongly affected by lateral variations in viscosity since the 1-D and 3-D predictions closely track one another (solid and dashed black line, Fig. 8e). We emphasize that the ice history (and relative size) of the Fennoscandian and Laurentide ice sheets, which are not explored here, will have a major affect on the GIA correction at these locations (Dendy et al. 2017; Rohling et al. 2017).

Ancient coral reefs in the Caribbean have long been used as palaeo sea level indicators. In particular, records from Xcaret on the Yucatan Peninsula (Blanchon et al. 2009) and various islands along the Bahamian archipelago (e.g. Hearty et al. 2007; Skrivanek et al. 2018; Dyer et al. 2021) have been influential due to the existence of fossil corals with low age uncertainty and good preservation. Being located on the tail end of the Laurentide peripheral bulge, these sites experience a small component of peripheral bulge subsidence (or equivalent sea level rise) in addition to continental levering. In both regions, the rate of sea level rise is higher in the 3-D simulation, which might be related to a low viscosity in the asthenosphere (see Section 3.2, Figs 6 and 5), a trend that is particularly notable for the Yucatan Peninsula (Figs 8c and d). Relative sea level predictions are slightly reduced in the 3-D GIA simulation that assumes a thinner lithosphere (blue dashed lines in Figs 8c and d). We also show predictions for Eleuthera in the northern Bahamas (yellow lines, Fig. 8d). For both 1-D and 3-D simulations, GIA predictions at Eleuthera are significantly different from the prediction at Great Inagua, which is expected given its location on the tail end of the peripheral bulge (Dyer et al. 2021). These differences demonstrate that applying a single GIA correction collectively to these sites is insufficient (Hearty et al. 2007; Clark et al. 2020) and that they should each be individually corrected prior to comparison (Dyer et al. 2021).

Far-field locations. The Seychelles and Western Australia are located in the far field of the former major ice sheets and have received substantial attention due to their high quality local sea level reconstructions (O’Leary et al. 2013; Dutton et al. 2015b). Our 1-D LIG sea level prediction in the Seychelles is relatively constant and slightly below the global mean. Incorporating lateral variations in viscosity leads to a slight upwards shift by 0.5–1.0 m, which is the result of a combination of a slightly thinner lithosphere and lateral viscosity perturbations (see Section 3.3; Figs 7 and 8a). The Western Australian coast is located on a hinge point, with higher sea level predicted offshore and lower sea level predicted on land when comparing 3-D and 1-D simulations (Section 3.3, Figs 5and 7). As a result, predictions using lateral variations in Earth structure are quite close and only slightly lower than predictions using 1-D Earth structure. At both far-field sites, relative sea level predictions are slightly increased at the beginning of the LIG and slightly decreased towards the end when assuming a 3-D earth model with a thinner lithosphere (blue dashed lines in Figs 8a and b).

4 DISCUSSION

The results presented above provide insight into the possible effects that lateral variations in lithospheric thickness and mantle viscosity can have on LIG sea level. Uncertainties remain in both the amplitude and pattern of viscosity perturbations, as we discuss in detail below. Ideally we would like to explore the full range of possible 3-D Earth structures, however, this is currently not computationally feasible. On the other hand, 1-D simulations are computationally inexpensive, and so we explore and discuss here two approximations: (1) we test whether the 3-D effects (i.e. the difference between a 3-D and 1-D simulation, where the spherical average of the former is given by the latter) are consistent for different choices of 1-D models and (2) whether 3-D GIA simulations can be approximated using 1-D simulations where the 1-D model differs from the spherical average of the 3-D earth model. We end our discussion by comparing our 3-D GIA predictions to relative sea level observations to understand how lateral variability in Earth structure may affect estimates of GMSL during the LIG.

4.1 What are the uncertainties in Earth structure?

There are three main factors that contribute to uncertainty in the mantle viscosity structure inferred from seismic tomography. The first involves the values of material properties that are used in the anelastic calibration (e.g. pressure- and temperature-dependence of the shear modulus and activation energy). The second is caused by intermodel differences in the seismic velocity structure imaged by different tomography studies. The third concerns the appropriate rheological deformation mechanism that is responsible for accommodating mantle flow during GIA. For the first two factors, our inverse anelastic calibration scheme provides a substantial advantage over traditional forward modelling approaches, which we illustrate in Fig. 9.

Figure 9.

Uncertainties in inferred Earth structure. (a) Standard deviation in inferred diffusion creep viscosity at 175 km depth for 1000 sets of anelastic parameters calibrated using the Richards et al. (2020) inversion scheme. Thick black line demarks lithosphere. (b) Same for a second suite of one thousand combinations of anelastic parameters, where each individual parameter is selected by randomly shuffling the values obtained in the construction of panel (a) and propagated into viscosity. Uncertainties are larger in this traditional forward mapping scheme due to the absence of information on the covariance between anelastic parameters. (c) Lateral viscosity perturbations at 175 km depth for an optimized calibration of CAM2016 seismic tomography model (Ho et al. 2016; Priestley et al. 2018). Thick black contour delineates the lithospheric portions. (d) As in (c), except for 3-D2015-07Sv model (Debayle et al. 2016). (e) As in (c), except for our preferred SL2013sv model (Schaeffer & Lebedev 2013) and its optimal value of activation energy, E_a = 304 kJ mol⁻¹. (f) Same as (e), except that lateral temperature variations obtained from the calibrated anelastic parametreization have been converted into viscosity using E_a = 560 kJ mol⁻¹, which is towards the upper end of the experimental range for dislocation creep in olivine.

The traditional approach is to adopt material properties that have been measured in laboratory experiments and convert, in a forward sense, from seismic velocity to temperature and viscosity. Including the inherent uncertainties associated with these measurements introduces a spread in inferred earth models. Our inverse calibration scheme, however, limits the number of acceptable combinations of material properties by retaining only those models that are consistent with the independent constraints on mantle structure (e.g. the thermal structure of oceanic lithosphere; Section 2.4.1). The approach reveals that there are trade-offs between the different material properties (Richards et al. 2020). Whilst uncertainty in any individual parameter remains large, exploiting their covariance results in a substantial reduction in the range of inferred earth models.

We illustrate this key benefit using a simple test. The initial step of the anelasticity optimization procedure is a coarse parameter sweep that is designed to locate the approximate position of the global misfit minimum. These parameters are then used as starting values in the second stage, which uses Powell’s algorithm to further minimize the misfit. For the test, we instead initiate this second stage from multiple different locations within the parameter space, discarding the result if the final misfit value returned by the algorithm is not smaller than the minimum value obtained in the coarse parameter sweep. In this manner, we obtain one thousand different sets of optimized anelasticity parameters that all yield satisfactory fits to the independent constraints. For each parameter, the range of optimal values across the one thousand sets is large and they remain individually uncertain. Nevertheless, the resulting standard deviation across all one thousand inferred viscosity structures is generally less than 0.2 orders of magnitude (Fig. 9a). Taking the set of one thousand values obtained for each individual anelastic parameter in the calibration stage, we can randomly shuffle them and construct a second suite of one thousand parameter combinations. This process yields the same total spread in individual material properties but removes information concerning their covariance (i.e. the information concerning which value of activation energy belongs with which value of reference viscosity, etc., is lost). Repeating the mapping from shear-wave velocity to viscosity (this time in a forward sense), we find that there is an approximately five-fold increase in the standard deviation of predicted viscosity models (Fig. 9b). Thus, our calibration scheme substantially reduces the uncertainty in inferred mantle viscosity structure. Exploiting parameter covariance in this manner is the strongest benefit of our inverse scheme over standard forward modelling practices.

The second source of uncertainty arises from differences in the starting seismic velocity structure between different tomography models. Choices including tomographic inversion technique, data content, reference velocity structure, and regularization all introduce inter-model differences. Traditional forward mapping schemes convert this variability into uncertainty in Earth structure. Our inverse calibration, however, reduces this uncertainty because it requires each tomography model to individually yield a temperature structure that is compatible with the independent constraints, thereby forcing some of the intermodel seismic velocity variation into the resulting optimal anelastic parameters. In Figs 9(c)–(e), we show results for three different surface wave tomography models, where the third case is the one used in this study: CAM2016 (Ho et al. 2016), 3-D2015-07Sv (Debayle et al. 2016), and SL2013sv (Schaeffer & Lebedev 2013). The resulting pattern of lateral viscosity perturbations is relatively consistent between the three models, although the features in SL2013sv tend to be slightly more localized and of higher amplitude in comparison to the other two. For each model, the root-mean-squared value of lateral viscosity perturbations outside of the lithosphere is 0.84, 0.75 and 0.78 Pa s, respectively. These values are more consistent with one another than the equivalent values obtained from forward mapping each tomography model using the same set of material properties (0.68, 0.60 and 0.78 Pa s).

The third source of uncertainty, that concerning the rheological mechanism by which the mantle deforms during GIA, is perhaps the most difficult to explore. The parametreization for anelasticity at seismic frequencies of Yamauchi & Takei (2016) yields a map of variations in the steady-state diffusion creep viscosity of the mantle. This deformation mechanism is consistent with our assumption in the GIA simulations that the mantle deforms like a Newtonian fluid. Nevertheless, it has also been suggested that deformation and flow during GIA may occur via dislocation creep, particularly in locations where strain rates are highest in the comparatively high-homologous-temperature asthenosphere (van der Wal et al. 2013; Huang et al. 2019). The dependence of dislocation creep viscosity on temperature (i.e. activation energy) has generally been found to be higher in laboratory experiments on olivine (430–570 kJ mol⁻¹ versus 240–425 kJ mol⁻¹ Karato & Wu 1993; Hirth & Kohlstedt 2003; Fei et al. 2012). Adopting mantle temperature variations obtained from the anelasticity parametreization, we see that applying this higher activation energy in the asthenosphere (between the lithosphere and 300 km depth) would lead to larger lateral variations in viscosity (Fig. 9f). We use this scenario to explore the impacts of larger lateral viscosity variations on our sea level reconstructions as might arise from dislocation creep, while still assuming diffusion creep in our calculations.

Fully propagating uncertainties in viscosity into our GIA predictions is outside the scope of this work. Nevertheless, in this and the following sections, we explore a few additional simulations. Using a viscosity model with a larger activation energy in the asthenosphere leads to submetre changes in predicted sea level in the far-field relative to a simulation with the defualt 3-D Earth structure (Fig. 10; dotted lines in Fig. 8). Note however that it can increase the overall 3-D effect, for example at the Seychelles, this simulation with a larger activation energy increases local sea level by approx. 1.5m relative to the 1-D simulation (black dotted and black solid line in Fig. 8a). In the near field, where sensitivity extends deeper into the mantle (where both models have the same viscosity variations), the difference in the two predictions is on the metre scale, which is smaller than the sea level change associated with introducing lateral variations in viscosity in the first place (Figs 5e and f).

Figure 10.

Effect of larger activation energy in the asthenosphere. A larger activation energy E_a leads to higher amplitude viscosity variations, which is in line with expectations for dislocation creep, although we do not explicitly model this rheology. (a and b) Difference in relative sea level between a 3-D viscosity model that uses a higher asthenospheric activation energy (above 300 km depth) and the reference case at the beginning (128 ka) and end (117 ka) of the Last Interglacial, respectively.

4.2 Are 3-D GIA effects dependent on the average viscosity profile?

Our default simulations (p55) use a lower mantle viscosity of 5 × 10²¹ Pa s, however Peltier et al. (2015) have argued for a weaker viscosity at this depth. To explore the dependence of 3-D GIA effects on the global average viscosity, we repeat our simulations with the VM5 viscosity profile in which lower mantle viscosity varies from 1.6 to 3.0 × 10²¹ Pa s (Peltier et al. 2015). GIA predictions using these two different 1-D viscosity profiles can differ significantly, especially in the near field (compare Figs 11a and b to Figs 4a and b). For example, the greater lower mantle viscosity in our default 1-D predictions (p55) results in lower sea level at Mallorca and higher sea level at Bermuda, Bahamas and Yucatan, compared to simulations using the VM5 viscosity profile (red lines in Fig. 8). Differences are smaller in the far field at sites such as the Seychelles and Western Australia (∼10s of centimetres).

Figure 11.

Effect of the choice of the 1-D background viscosity on LIG sea level. (a and b) Relative sea level at the beginning (128 ka) and end (117 ka) of the Last Interglacial, assuming the radially symmetric Earth structure VM5 (Peltier et al. 2015). (c and d) Same as (a) and (b) but including lateral variations in lithospheric thickness and mantle viscosity. (e and f) Differences in relative sea level between model simulations that do and do not account for lateral variability in Earth structure. (g and h) A comparison of the effect of lateral variations in viscosity when superimposed on the VM5 and p55 viscosity profiles. Plots show the difference in the 3-D effect, that is the difference between panels (e) and (f) of this figure and panels (e) and (f) of Fig. 5.

Using these results, we investigate whether the incorporation of lateral variations in viscosity has the same effect whether the p55 or VM5 depth average viscosity is adopted in the simulation. Note that the lateral variations in viscosity shown in Fig. 3 are superimposed on these two 1-D profiles such that the spherical average of the logarithm of viscosity at each depth remains unaffected. Comparing the results in Figs 5(e) and (f) to Figs 11(e) and (f) indicates that the impact of lateral viscosity structure is qualitatively similar; however, differences in magnitude and geometry exist (Figs 11g and h). This similarity is also evident when comparing results at specific locations: At Mallorca, for example, we find that while the choice of the 1-D profile results in two different sea level predictions (black versus red solid line in Fig. 8e), the signal due to the introduction of lateral variations is consistent (black and red dashed lines). On the peripheral bulge of the former Laurentide ice sheet, this signal has the same sign but differs in magnitude from site to site and is generally larger when adopting the VM5 viscosity profile. For example, at Bermuda, the sea level predictions based on the p55 1-D and 3-D simulations are similar (within ∼0.5 m; black solid versus dashed line in Fig. 8f), while the effect of adding lateral variations in viscosity is much larger when assuming the VM5 viscosity profile (+3 m towards the end of the LIG; red solid versus dashed line in Fig. 8f). For the three sites in the vicinity of the Laurentide peripheral bulge (Xcaret, Bahamas and Bermuda), the two 3-D predictions are more consistent with one another than their associated 1-D predictions, particularly for Xcaret where these differences remain less than 1 m throughout the LIG (Figs 8c, d and f). In the far field, we find that the introduction of lateral variations in viscosity tends to consistently increase local relative sea level predictions at the Seychelles by up to ∼1 m and decrease them by a similar amount in Western Australia.

4.3 Can 3-D simulations be approximated with 1-D simulations that are not the spherical average of the 3-D earth model?

Given the computational expense of 3-D GIA simulations, it is worth investigating whether a simulation with a suitable 1-D viscosity profile, which is not necessarily the spherical average of the 3-D earth model, can be used to approximate the 3-D result with sufficient accuracy. Powell et al. (2019) considered synthetic GPS observations in Antarctica and found that 1-D simulations tuned to Earth structure local to the sites do not provide consistently accurate approximations to the 3-D synthetic predictions. Hartmann et al. (2020) have proposed an approach which combines the result from different 1-D simulations to approximate the 3-D result. They focus on Antarctica and argue that the approach has promise, but concede that it might be inaccurate in areas where the ice load and sea level observation are relatively distant from one another. The latter situation is the case for most sea level studies that consider observations distant from the former ice margins, such as this study. Crawford et al. (2018) used an adjoint approach to produce 3-D sensitivity kernels that isolate regions of the mantle that are sampled by a given sea level record and predict whether increasing or decreasing viscosity in these regions will lead to a better fit between the model prediction and the sea level observation. They found that the sensitivity is centred below the location of the sea level record and extends towards the locations of ice melt. Moreover, their time-dependent sensitivity kernels indicate that the region of greatest sensitivity will vary over time. These results suggest that approximating 3-D Earth structure using 1-D simulations may be challenging. Nevertheless, we explore two approaches here: First, we take a depth-average of the 3-D earth model in the vicinity of each individual sea level site (averaging is performed in logarithmic space and within a maximum distance of 3° around the sea level site; Figs 12a–f) and repeat our 1-D simulations using this local structure. Secondly, we use a broad suite of different viscosity profiles to assess whether any of them can provide a good approximation to the 3-D GIA result.

Figure 12.

Approximating 3-D Earth structure with radially symmetric structure. (a–f) Locally averaged viscosity structure as a function of depth from the 3-D earth model (purple), global average viscosity (black), and 3-layer 1-D viscosity profile that best fits the 3-D GIA prediction at each location (green). Yellow band in panel (a) shows the full range of 1-D models explored here. To compute the locally averaged viscosity structure, the viscosity below the site was averaged across a maximum distance of 3° from the sea level site at each depth. (g–l) Local relative sea level predictions at selected sites (see caption of Fig. 8 for exact locations) using the 1-D viscosity profiles shown in panels (a)–(f) with the same colours and also including predictions for the full 3-D Earth structure (black dashed line). Thick grey line is the global mean sea level value assumed in the GIA model. Light green range marks the 1σ uncertainty range for the ensemble of 1-D runs explored here. (m–r) Parameter sweeps through upper and lower mantle viscosity (see text) at optimal lithospheric thickness for each site, showing misfit between each individual 1-D prediction and the 3-D prediction. The optimal lithospheric thickness is noted in the bottom right corner of each panel. The earth model with the minimum misfit is shown by the white circle (this model is given by the green line in panels (a)–(f).

The locally averaged Earth structure obtained from the 3-D model below the six key sea level locations shows that most of them have a weaker than average upper mantle viscosity (Figs 12a–f), which is not surprising given that most of them are distant to subduction zones and cratonic regions. The only exception is Mallorca, which has a higher than average viscosity in the deeper half of the upper mantle due the presence of a subducted slab. Local viscosity variations in the lower mantle are more variable, with larger differences in particular in the vicinity of the core-mantle boundary, which will have limited influence on the GIA response. Figs 12(g)–(l) compare the result from the 1-D simulations adopting local Earth structure (purple line) with the full 3-D result (black dashed line). The two results are consistent for the Yucatan peninsula and Western Australia (Figs 12g and i), but do not agree well elsewhere.

We next test a range of 1-D earth models to investigate which (if any) structure approximates the local 3-D result for each site. We consider 48 different three-layer radial earth models that each consist of an elastic lithosphere overlying isoviscous upper and lower mantle regions. We systematically vary upper and lower mantle viscosity across 3–5 × 10²⁰ Pa s and 3–40 × 10 ²¹ Pa s, respectively, and test two different lithospheric thicknesses (71 and 96 km; see Fig. 12a). The 1σ range of all model simulations is shown in green in Figs 12(g)–(l). Sites in the far-field are most sensitive to lithospheric thickness variations and upper mantle viscosity since continental levering is an important driver of sea level change for these sites. Sites on the peripheral bulge of the former ice sheet are more sensitive to mantle viscosity: Mallorca is most sensitive to lower mantle viscosity, and Bermuda, Great Inagua and Xcaret are equally sensitive to upper and lower mantle viscosity.

We next compare our predictions for each 1-D simulation to the 3-D result at the six sites, calculating misfit using the root-mean-square difference in relative sea level over the LIG (between 117 ka and 128 ka). The misfit, which is shown as a function of upper and lower mantle viscosity in Figs 12(m)–(r), shows a strong dependence on Earth structure for near-field sites (Figs 12o–r) and a weaker dependency for far-field sites (Figs 12m and n). We find that the best fitting 1-D earth model at each site does produce a sea level prediction that matches the 3-D simulation reasonably well (green line compared to black dashed line in Figs 12g–l). It is difficult to compare the local Earth structure to the best-fitting 1-D Earth structure given the coarse resolution of the latter, however, the two show some consistency at far-field sites (green line compared to purple line in Figs 12c–f). Differences between the local and best fitting 1-D Earth structure are expected given the broad sensitivity of sea level observations, which integrates Earth structure across wide regions of the mantle (Crawford et al. 2018). While the difference in the relative sea level prediction using a local structure versus the full 3-D Earth structure argues against using the former as an approximation for the 3-D result, the suite of 1-D results suggest that a suitable 1-D approximation may exist at each site. Inferring radial Earth structure from observations at these sites would generally lead to overestimates of the average lower mantle viscosity. In other words, while the global average viscosity in the lower mantle in these simulations is 5 × 10²¹ Pa s, the viscosity of the best fitting 1-D model at all sites except for the Seychelles is higher than that value.

Finally, the above analysis raises the question: how useful is it to use a range of 1-D viscosity models when estimating uncertainties in the GIA correction (particularly uncertainties introduced by Earth structure)? The green band in Figs 12(g)–(l) shows the 1σ uncertainty range associated with the full ensemble of 1-D earth models used here (the mean is not shown, but it sits in the middle of the light green band). In the near field, the 3-D result falls within the 1 σ range of 1-D predictions (Figs 12i–l). In the far field, the uncertainty range is relatively narrow and the 3-D prediction falls just outside of this 1σ range, but within the 2σ range (Figs 12g and h, note that the 2σ is not shown). We thus consider that results based on a range of 1-D model runs may provide a suitable estimate of the uncertainty associated with the potential signal from lateral variations in viscosity structure.

4.4 How do lateral variations in Earth structure affect estimates of LIG global mean sea level?

Estimates of GMSL during the LIG are based on sea level observations (such as corals or speleothems) from this time period. The locations we have chosen for our investigation (Figs 8 and 13) are among the sites with the most reliable local sea level records. The inferred GMSL estimate at each site is given by the difference between observed sea level and that predicted by the GIA simulation. Sea level during the LIG will also vary spatially depending on which ice sheet is driving the excess melting (Hay et al. 2014), an issue which is not explored here.

Figure 13.

GIA time series at key sites with LIG sea level records. (a-f) Relative sea level at the Seychelles (La Digue), Western Australia (Cape Range), Yucatan Peninsula (Xcaret), Bahamas (Great Inagua and Eleuthera, Whale Point), Mallorca (Coves del Pirata), and Bermuda (Grape Bay), respectively. Thick grey line is the eustatic (global mean sea level) value assumed in the GIA model. Thus, any predicted relative sea level change during the LIG is only due to GIA and not global mean sea level changes. Black lines show predictions of relative sea level using a 1-D (solid line) and the 3-D (dashed line) GIA model with spherical average given by the p55 viscosity profile. Grey lines are results for different 1-D and 3-D earth models from Fig. 8. Light blue lines show the inferred relative sea level at each site based on a variety of observations, with shaded regions marking uncertainties cited in the original publications: Seychelles (Dutton et al. 2015b); Western Australia, line without uncertainty (O’Leary et al. 2013) and line with uncertainty (Dutton & Lambeck 2012); Yucatan (Blanchon et al. 2009); Bahamas, with 1σ and 2σ uncertainty (Dyer et al. 2021); Mallorca (Polyak et al. 2018); and Bermuda, time-varying prediction (Hearty 2002) and constant prediction based on the highest reported Devonshire marine member which has large age uncertainties (Muhs et al. 2020). (g) Locations from (a)-(f) shown on a map.

In Mallorca, phreatic overgrowths on speleothems (POS) have been used to reconstruct local sea level, which the authors infer to be relatively stable throughout the LIG (Polyak et al. 2018, Fig. 13e). Given that relative sea level is predicted to steadily rise due to GIA, Polyak et al. (2018) concluded that GMSL must be falling over the LIG in order to result in constant relative sea level. Our result indicates that an even greater fall of sea level would be required if lateral variations in viscosity are accounted for (Fig. 13e). Bermuda is the other near-field site in our analysis, and stratigraphic and coral evidence suggests that local sea level peaked around 6–8 m above present; however the exact timing and evolution is controversial due to insufficient age control (Hearty 2002; Muhs et al. 2020). Accounting for lateral variations in viscosity will tend to reduce the magnitude of the inferred GMSL and, assuming that the highstand was recorded late in the LIG, implies only a few metres of excess GMSL during that time.

More distal near-field records from the Yucatan Peninsula and the Bahamas show locally rising sea level, which are recorded by extensive coral reefs. At the ecological park of Xcaret, Blanchon et al. (2009) identified a lower and upper reef crest (Fig. 13c). In the Bahamas, Dyer et al. (2021) used coral and sedimentary evidence combined with a large suite of radially symmetric GIA models to calculate a posterior relative sea level history that exhibits an early sea level rise, followed by slightly falling sea level before culminating in a final rise (Fig. 13d). This history is in agreement with earlier analyses from this location (Dutton & Lambeck 2012; Skrivanek et al. 2018). If one were to assume 1-D Earth structure at these locations, one would infer ∼ 3–4 m of excess GMSL early and a smaller excess late in the interglacial, with a GMSL lowstand in the interim. The 3-D GIA predictions are higher than the 1-D predictions towards the end of the LIG, which may lower the inferred GMSL at the end of the LIG.

Far-field records along the western coast of Australia and in the Seychelles are also based on coral outcrops. O’Leary et al. (2013) dated corals at several locations in Western Australia and inferred an early rise in local sea level that was followed by a GIA-driven sea level fall, which resulted in erosion of a coral platform (Fig. 13b). Additional higher corals were interpreted to reflect a late rise in sea level. Dutton & Lambeck (2012) inferred a similar planated surface during the first half of the LIG and interpreted the higher corals to be tectonically deformed (see also Sandstrom et al. 2020). GIA and therefore inferred GMSL would be marginally impacted by 3-D Earth structure at this location, which tends to increase inferred GMSL (by ∼ 0.5 m). Inferred GMSL is 3–4 m at the beginning of the LIG (in line with earlier estimates) and remains at that level to the end of the LIG if the high corals are discounted or increases to ∼9 m if they are not. Finally, extensive coral reefs are absent on the Seychelles, but individual corals and coralline algae are attached to granitic bedrock (Dutton et al. 2015b; Vyverberg et al. 2018) and found at high elevations, leading to an interpreted local sea level of around 6–7 m above present early in the LIG (Dutton et al. 2015b). 3-D GIA results tend to increase the predicted relative sea level, which decreases the inferred GMSL (Fig. 13a). The magnitude of this effect ranges from 0–1.5 m (where this range includes the simulation with a larger activation energy in the asthenosphere) leading to an inferred GMSL early in the LIG that remains larger than that at most other sites (6–8.5 m). Increasing the GIA prediction for local relative sea level (and hence reducing the inferred GMSL) is possible by decreasing the lithospheric thickness in this region (Fig. 7), which enhances continental levering. However, this effect would also be expected to occur during the Holocene and could result in an early Holocene sea level highstand (depending on the GMSL history), which has not been observed (Woodroffe et al. 2015).

Inferences of GMSL during the LIG described above are based on a limited number of 3-D simulations, and a rigorous analysis would require testing a significantly larger suite of Earth structures. In addition to Earth structure, there are several major uncertainties associated with the ice history that are not explored in this analysis, but will be briefly summarized: (1) the calculations performed here begin at 150 ka (Fig. 1), which assumes that the ice-Earth system was in isostatic equilibrium at this time. We have performed 1-D simulations that include earlier glacial cycles and found that this effect can cause differences on the order of 1 m in areas of the peripheral bulge and smaller (decimetre scale) in the far field; (2) GIA across the LIG will be sensitive to the specific ice sheet configuration adopted during the penultimate deglaciation, an uncertainty explored in detail elsewhere (Dendy et al. 2017; Rohling et al. 2017), and this factor will be particularly crucial to consider when attempting to reconcile relatively near-field sites such as Mallorca and Bermuda; (3) GIA predictions of relative sea level during the LIG are also sensitive to the ice history during the last glacial cycle (Lambeck et al. 2012). Here, we have assumed that sea level was relatively high during MIS 3 due to a small Laurentide ice sheet, following the results of Pico et al. (2017). If we were to assume that the Laurentide ice sheet was larger during MIS 3, it would lead to a further increase in predicted relative sea level during the LIG at sites close to the former Laurentide ice sheet; and (4) Ice melt during the LIG will drive spatially variable sea level changes, and this should be accounted for when comparing GMSL inferences from different locations (Hay et al. 2014).

5 CONCLUSION

In this study we describe GIA predictions based on a new model of Earth’s 3-D viscoelastic structure inferred from recent global tomographic models (Schaeffer & Lebedev 2013, 2014; French & Romanowicz 2014). We use an upper mantle anelastic parametreization that relates shear wave speed to diffusion creep viscosity and is based on laboratory deformation experiments (Yamauchi & Takei 2016). The parameters within these relationships are calibrated such that the resulting temperature variations match a series of independent observables (Richards et al. 2020) and this reduces the uncertainty in the inferred viscosity. We note that the apparent viscosity over ice age timescales might deviate from the steady-state viscosity due to transient behaviour (Lau & Holtzman 2019), which is not explored here.

We use this new model of Earth’s internal structure to produce the first estimates of GIA-driven sea level change across the LIG that incorporate lateral variations in viscoelastic structure. We find that GIA predictions of relative sea level based on 3-D versus 1-D Earth structure have metre-scale differences in both the near and far field. We explore the mechanisms responsible for these differences and demonstrate how effects such as forebulge dynamics and continental levering are influenced by the presence of lateral variations in lithospheric thickness and underlying mantle viscosity. A more detailed examination of these differences is possible using 3-D sensitivity kernels (Al-Attar & Tromp 2014; Crawford et al. 2018).

The effect that lateral viscosity variations have on sea level is weakly dependent on the globally averaged 1-D viscosity structure that these variations are superimposed on: Using two different 1-D profiles, we find that the difference between 3-D and 1-D predictions of LIG sea level differ more in magnitude than in geographic pattern. Thus, our results can be cautiously used as a first-order guide to whether lateral mantle viscosity variations might increase or decrease relative sea level in comparison to 1-D GIA predictions.

Given the computational expense of 3-D GIA simulations, it is important to consider if and how well such simulations can be approximated by 1-D GIA modelling. We find that 1-D simulations that assume local Earth structure within a 3° radius of the site do not produce results that are representative of the 3-D result, which is consistent with earlier findings (Powell et al. 2019; Hartmann et al. 2020). However, a suite of 1-D simulations suggests that a suitable and unique 1-D approximation may exist for each site and we speculate that appropriate values for such a model might be found by averaging 3-D structure over mantle regions characterized by high sensitivity (Crawford et al. 2018).

Lastly, we compare our predictions of GIA for 3-D earth models to local LIG sea level reconstructions to investigate the implications of such models for estimates of GMSL during the LIG. It is noteworthy that lateral variations in mantle viscosity perturb predictions in a manner that may help to reconcile the mismatch in inferred GMSL early in the LIG from the Seychelles, where they lower this value, and Western Australia, where they increase it. However, this effect is not large enough to fully bring published estimates from these two sites into accord. Our results show that lateral variations in Earth structure are important to consider when reconstructing past sea level and estimating peak GMSL (or minimum ice volumes) during periods of relative ice age warmth.

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

APPENDIX A: ANELASTIC PARAMETREIZATION

The anelastic parametreization of Yamauchi & Takei (2016) represents linear viscoelasticity in the frequency domain using a complex compliance, the real component of which refers to the amplitude of strain that occurs in phase with the driving stress and the imaginary component refers to strain that is π/2 radians out of phase (and gives rise to dissipation). Processes responsible for accommodating anelastic deformation are represented using a relaxation spectrum that consists of a high-frequency peak superimposed on top of a monotonic background. The height and width of the high-frequency peak is a function of the material’s homologous temperature, |$\frac{T}{T_m}$|⁠. The background takes advantage of the Maxwell frequency ‘master variable’ scaling results of McCarthy et al. (2011), which crucially showed that the effects of variations in composition, grain size and temperature on attenuation can be accurately predicted using the corresponding effect of these parameters on the steady-state diffusion creep viscosity of the material.

In their laboratory experiments on organic rock-analogues, Yamauchi & Takei (2016) were able to ascertain the values of several anelastic parameters that are thought to be consistent between different polyscrystalline materials (e.g. the centre frequency of the relaxation peak). Seven other parameters are specific to each individual material and need to be independently determined. These include the unrelaxed shear modulus at reference conditions, its dependence on temperature and pressure, the viscosity at reference conditions (1200 °C and 1.5 GPa), its dependence on temperature and pressure (activation energy and activation volume, respectively), and the solidus gradient. Assuming that suitable values for these parameters can be ascertained, the parametreization allows shear wave velocities to be converted into temperature and steady-state diffusion creep viscosity in a self-consistent manner. The exact form used in this study is given by eqs (3)–(17) of Richards et al. (2020).

ACKNOWLEDGEMENTS

This work has been supported by the National Science Foundation under grant OCE 18-41888 and ICER 19-28146 (JA) and OCE 17-02684 (JXM). MH and JXM acknowledge support from the National Aeronautics and Space Administration (grant NNX17AE17G). MH acknowledges support from the Australian government’s Exploring for the Future program. FDR acknowledges support from the Schmidt Science Fellows program, in partnership with the Rhodes Trust, and the Imperial College Research Fellowship scheme. We are grateful to Tamara Pico for providing her ice reconstruction for MIS 3 to MIS 5d, Sia Ghelichkhan for discussions concerning conversion of seismic velocities into Earth structure, and Andrew Lloyd for discussions about sensitivity kernels and general feedback. The authors acknowledge PALSEA, a working group of the International Union for Quaternary Sciences (INQUA) and Past Global Changes (PAGES), which in turn received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences.

AUTHOR CONTRIBUTIONS

JA designed the research and performed the sea level analysis with input from all authors. MH and FDR performed conversions from shear wave speed to temperature and explored related uncertainties. KL performed the 3-D GIA calculations and JA performed the 1-D GIA calculations. JA and MH wrote the majority of the manuscript with input from all authors.

DATA AVAILABILITY

The data underlying this paper are available as supporting information with the online version of this manuscript and in the accompanying OSF database (https://osf.io/m9qcz/).