Glacial hydro-isostatic adjustment at mid-ocean ridges

SUMMARY

Recent studies have suggested a link between ice age sea level fluctuations and variations in magma production and crustal faulting along mid-ocean ridges based on the detection of Milankovitch cycle frequencies in topography off several ridges. These fluctuations have also been connected to variability in hydrothermal metal fluxes near ridges. Ice age sea level calculations have shown that the sea level change across glacial cycles will be characterized by significant geographic variability, that is, departures from eustasy, due to the gravitational, deformation and rotational effects of the glacial isostatic adjustment (GIA) process. Using a state-of-the-art GIA simulation that incorporates 3-D variations in Earth viscoelastic structure, including plate boundaries and updated constraints on the magnitude and geometry of ice mass fluctuations, we predict global sea level changes from last glacial maximum (LGM, 26 ka) to present and from the penultimate glacial maximum (143 ka) to the last interglacial (128 ka). We focus on the results along three ridges: the Mid-Atlantic Ridge, Juan de Fuca Ridge and East Pacific Rise, which are examples of slow, intermediate and fast spreading ridges, respectively. Sea level change across the Mid-Atlantic Ridge shows the greatest variability, ranging from a sea level fall greater than 200 m in Iceland to a maximum rise of ∼150 m in the South Atlantic, with significant non-monotonicity north of the Equator as the ridge weaves across the field of sea level changes. We also calculate changes in crustal normal stress from LGM to present-day across the Mid-Atlantic and Juan de Fuca Ridges and the East Pacific Rise. These results indicate that the contribution from ice mass changes to the crustal stress field can be significant well away from the location of ancient ice complexes. We conclude that any exploration of the hypothesized links to magma production and crustal faulting must consider both ocean and ice loading effects and, more generally, the profound geographic variability of the GIA process.

1 INTRODUCTION

Mass transfer between grounded ice sheets and glaciers and the global ocean during the Late Pleistocene ice age has been connected to the frequency of volcanic eruptions, whether as a response to changes in crustal stress due to ice unloading on land in the near field of glaciation (e.g. Sigvaldason et al. 1992; Jull & McKenzie 1996; Jellinek et al. 2004; Praetorius et al. 2016; Rawson et al. 2016) or to sea level changes in the far field of the cryosphere (e.g. McGuire et al. 1997; Kutterolf et al. 2013; Schindlbeck et al. 2018). In these studies, the increase in volcanism was found to follow, with relatively minor time lag, the period of greatest reduction in normal stress (e.g. Kutterolf et al. 2013), and tephra records which extend over multiple glacial cycles were shown to have statistically significant Milankovitch periodicities (e.g. McGuire et al. 1997; Kutterolf et al. 2013; Schindlbeck et al. 2018). A variety of physical mechanisms have been proposed that connect crustal stress changes and magmatic output. Studies of Iceland, for example, argued that the reduction in the crustal stress driven by ice unloading leads to either a greater sampling of the magma chamber (Gudmundsson 1986) or to greater magma production (Jull & McKenzie 1996).

The connection between volcanism and ice age loading has also been suggested for mid-ocean ridges (MORs), where rising sea level (deglaciation) may hinder volcanic output, particularly at slow spreading ridges, and falling sea level (glaciation) may amplify this output (Huybers & Langmuir 2009). This connection was initially supported by the detection of Milankovitch periodicities in abyssal plain topography (Huybers & Langmuir 2009), but subsequent efforts to detect this topography signature have yielded mixed results (Crowley et al. 2015; Goff 2015, 2020; Tolstoy 2015; Boulahanis et al. 2020). An analysis by Lund and Asimow (2011) reinforced the argument that sea-level variations can drive observable fluctuations in abyssal plain topography, and they suggested other potential proxies for magma-flux, including oceanic radiocarbon, the oceanic Osmium isotopic ratio, and hydrothermal particle flux; the latter has been detected in metal deposits across one or more glacial cycles at sites in the northern Mid-Atlantic Ridge (Middleton et al. 2016), the Juan de Fuca Ridge (Costa et al. 2017) and the East Pacific Rise (Lund et al. 2016, 2019). These studies point to possible feedback between sea level change and CO₂ emissions that may act as a driver for ice age climate (Huybers & Langmuir 2017). However, other analyses have argued that variations in magma production would be insufficient to produce significant topographic variability (Olive et al. 2015; but see Huybers et al. 2016) and, further, that the characteristic width of abyssal hills as a function of spreading rate is consistent, instead, with fault-generation in response to sea level changes (Goff 2015; Olive et al. 2015). The latter was supported by Huybers et al. (2022), though they also argued that at intermediate spreading rates the trend in abyssal hill widths supported the connection between sea level changes and magma production. They also describe a faulting model consistent with these nuanced observations.

The above studies suggest—albeit with significant, ongoing challenge and debate (Goff 2015, 2020)—that changes in sea level across the Late Pleistocene glacial cycles may imprint on magma production, hydrothermal fluxes and fault generation at MORs and the adjacent abyssal hills and, more generally, may drive variability in ice age climate through CO₂ outgassing. However, two important issues have been neglected in studies supporting these connections. The first is that sea level changes driven by ice age cycles are characterized by significant geographic variability that reflect the deformational, gravitational and rotational effects of the glacial isostatic adjustment (GIA) process (e.g. Milne & Mitrovica 2008; Lambeck et al. 2014). Second, variations in normal stresses at MORs may also be impacted by isostatic adjustment driven by ice mass changes since this deformation field can extend at great distance from the location of the major Late Pleistocene ice complexes (Milne & Mitrovica 2008; Coulson et al. 2021). In this paper, we use a recent, state of the art model of the GIA process to explore how these sources of variability are manifest along the Earth's networks of MORs, with a particular focus on the fast-spreading East Pacific Rise, the slow-spreading Mid-Atlantic Ridge and the intermediate-spreading Juan de Fuca Ridge. We begin by demonstrating the ice-age cycle variability in sea level change along the Earth's networks of MORs. Next, we compute near-surface normal stresses at a subset of these ridges and partition this stress field into contributions from ice and ocean loading. At some level, variations in GIA-induced crustal stresses—whether driven by meltwater or ice loading—would modulate the suite of hypothesized connections described above.

2 METHODS

We compute gravitationally self-consistent sea level and crustal stress variations on a deforming and rotating viscoelastic Earth model following the approach described in Kendall et al. (2005). The modelling considers shoreline migration, water redistribution in areas where grounded, marine based ice sheets advance or retreat, and the feedback into sea level of rotation axis reorientation. The GIA simulations are performed using the SeaKon finite volume software described in Latychev et al. (2005). The spatial resolution varies from 10 km near the surface of the model to ∼50 km close to the core–mantle boundary. Two inputs are required in the simulations, an Earth model and global ice history, and we describe each in turn.

The vast majority of previous GIA studies have adopted 1-D viscoelastic Earth models in which viscosity structure varies only with depth. In contrast, we adopt a new, 3-D viscoelastic model constructed by Austermann et al. (2021) from the work of Richards et al. (2020). The model is consistent with a range of results from seismic tomography and attenuation measurements, high-pressure laboratory results on mantle materials and geological data. The density and elastic structures of the model are taken from the seismic model PREM (Dziewonski & Anderson 1981).

We adopt the PC2T ice history developed over the last decade in a series of articles and based on a wide range of geological markers of sea level change (Creveling et al. 2017; Dendy et al. 2017; Pico et al. 2017; Austermann et al., 2021). The model extends back to 160 kyr before present (i.e. 160 ka) and includes two full deglaciation phases with distinct geometries of ice melt, both characterized by a global mean sea level (GMSL) rise of ∼130 m (Fig. 1). The most recent of the two, extending from Last Glacial Maximum (LGM, 26 ka) to ∼5 ka, is largely based on the ICE-6G model of Peltier et al. (2015), and is characterized by ice cover over North America (including the Laurentide, Cordilleran and Innuitian Ice Sheets) that accounts for ∼80 per cent of the GMSL rise. The earlier deglaciation, extending from the penultimate glacial maximum (PGM, 143 ka) to the last interglacial (LIG, 128 ka), is a more rapid event. Moreover, the model adopts the PGM ice geometry in the Northern Hemisphere of Lambeck et al. (2006) in which the ice sheets over North America plus Greenland and over Eurasia contribute a nearly equal amount (∼56 m) to the subsequent GMSL rise. (A side-by-side comparison of Northern Hemisphere ice cover during the PGM and LGM is provided by Dendy et al. 2017; Fig. 3.) Computing the global sea level change across these two deglaciation phases thus provides a measure of the sensitivity of our sea level predictions at MORs to the geometry of ice mass flux from times of glacial maxima to subsequent interglacials.

3 RESULTS

3.1 Sea level changes—a global view

The predicted change in sea level from LGM to present day computed with the PC2T ice history and the 3-D viscoelastic Earth model is shown in Fig. 2. The GMSL rise across this period is ∼125 m (Fig. 1). The basic physics underlying LGM to present-day sea-level change was discussed by Milne & Mitrovica (2008). In the near field of the location of the ancient ice sheets, the predicted sea level fall can reach many hundreds of meters and is well off the scale of the plot. This sea level fall is largely due to two effects: the migration of water away from the diminishing ice sheets as they lose gravitational pull on the water and viscoelastic post-glacial rebound of the crust. This near field region extends somewhat beyond the maximum perimeter of the ice sheets. At greater distance, in the so-called far field, the magnitude of sea level rise increases to a maximum of ∼150 m, well above the GMSL rise, in the vicinity of South America.

Secondary processes impacting sea level are also visible on the figure. Close to far-field shorelines, the predicted sea level rise decreases moving toward the continent in response to crustal uplift due to ocean loading in a process called continental levering (Nakada & Lambeck 1989). The localized zone of higher sea level rise (deep blue) to the northeast and directly south of Iceland, as well as the more muted rise southwest of Iceland, is due to the combined subsidence of the crustal bulges at the periphery of the Laurentide, Fennoscandian and Greenland Ice Sheets. The same localized zone of higher sea level rise is also evident off the west coast of Canada, at the periphery of the Cordilleran Ice Sheet. Peripheral bulges encircled all glaciation centres, but the localized sea level rise associated with this subsidence is often masked by the migration of water away from the near field due to gravitational effects. Finally rotational feedback on sea level is evident in the skewing of the regions of maximum sea level rise (e.g. from the southeast Pacific to the northern Indian Ocean).

3.2 Sea level changes—mid ocean ridges

Fig. 3 highlights results for the prediction in Fig. 2 along the Mid-Atlantic Ridge (left frame, and right frame—blue points). As one moves north along the ridge, the prediction reaches a maximum of ∼150 m (point A) and falls relatively monotonically until a latitude of ∼30° N (point B). As the ridge veers eastward, a local minimum occurs at point C. Further north along the ridge, the prediction shows relatively large oscillations. A dip in the sea level prediction occurs as the ridge samples the near field, and gravitational effects, of the Laurentide and Greenland Ice Sheets (point D). Next a local maximum occurs on the subsiding bulges at the periphery of various zones of deglaciation (point E), and then the prediction incurs a rapid drop toward Iceland as the gravitationally driven migration of water away from the melt zone and crustal uplift due to unloading begin to dominate. Point F lies just south of Iceland and the net sea level change from LGM to present day at the site is +15 m. Within Iceland, sea level is predicted to fall by up to ∼260 m. North of Iceland the ridge samples the bulge adjacent to ice melting over Iceland, Greenland and Fennoscandia, and the associated crustal subsidence dominates other effects to produce a local maximum in the prediction (point G). Finally, as the ridge veers toward Fennoscandia and Svalbard, gravitational effects dominate, and the predicted sea level change falls just below +30 m (point H), then recovers to ∼100 m as the ridge veers away from this near field signal (point I). In summary, the prediction shows large variations from the GMSL rise in the PC2T ice history of ∼125 m (Fig. 1) associated with the net ice mass flux from LGM to present day.

Fig. 4 highlights sea level change predictions along the East Pacific Rise, which are characterized by a relatively monotonic fall in magnitude as one moves along the ridge from a site in the southeast Pacific (point A, ∼150 m) to the southern end of Baja California (point B, ∼100 m). As is clear from Fig. 2, this trend largely reflects the pattern of water migration away from the major centres of glaciation in the Northern Hemisphere. Once again, the predicted sea level change shows large departures from the GMSL rise of 125 m across this time interval.

Analogous results for the Juan de Fuca ridge are shown in Fig. 5. Here, the variability in the sea level signal is more muted, as the ridge lies in a zone in which crustal subsidence at the periphery of the Cordilleran Ice Sheet is largely compensated by the gravitational migration of water away from the ice sheet (Figs 2 and 5). As one moves north along the ridge, from point A to B, there is a general trend downward in the predicted sea level change, from 110 to 93 m, as one samples sites closer to areas in which the gravitational effect begins to dominate (left, Fig. 5, green–yellow–red contours).

Figs S2 and S3 (Supporting Information) explore the sensitivity of the predictions to variations in Earth structure using a suite of 1-D Earth models. Fig S2 (Supporting Information) reproduces the LGM-to-present results in Figs 3–5 and compares these to calculations based on a 1-D Earth model identical the spherically averaged profile of the 3-D model. The 1-D results are significantly smoother than the profiles based on the 3-D Earth model but exhibit similar geographic trends at the largest spatial scale. Fig. S3 (Supporting Information) considers a range of 1-D Earth models in which the lithospheric thickness, upper mantle viscosity and lower mantle viscosity are varied in the GIA simulation. The sea level predictions across the Mid-Atlantic Ridge south of Iceland and the East Pacific Rise are most sensitive to lower mantle viscosity, reflecting, in part, the broad spatial extent of the ocean and Laurentide ice loading. As this viscosity is increased in the modelling the magnitude of the geographic trend in these predictions increases. North of Iceland, the predictions for the Mid-Atlantic Ridge show greatest sensitivity to the adopted lithospheric thickness. Finally, results for the Juan de Fuca Ridge are sensitive to all aspects of the Earth model, and most notably the lithospheric thickness. These sensitivities are a consequence of the proximity of the ridge to the Cordilleran Ice Sheet.

Fig. 6 shows a representative sample of relative sea level (RSL) curves along each of the three ridges sampled in Figs 3–5. The curves for sites on the East Pacific Rise mimic the GMSL curve for the PC2T ice history (Fig. 1) with amplitudes monotonically increasing as one moves south along the ridge. The RSL curves along the Juan de Fuca Ridge show a similar time history, although with more variability, particularly during the LGM and PGM. This latter complexity reflects the proximity of sites D and F to the ancient Cordilleran Ice Sheet. Finally, the RSL curves for sites along the Mid-Atlantic Ridge show major departures from the GMSL curve that increase as one considers sites closer to Iceland. The RSL curve at Site F, just south of Iceland, shows little resemblance to the GMSL curve of Fig. 1, with a variation that lies within 20 m of present-day level across the entire time history of the simulation and periods of sea level fall during the two deglacial phases (e.g. at 15–7 and 135–128 ka)

As noted in the Methods section, the geometry of ice mass flux across the penultimate deglaciation phase (143–128 ka) was notably distinct from the geometry in the post-LGM phase, with the former characterized by comparable reductions in ice volume over North America and Eurasia. The PC2T ice history covers both periods and, to facilitate comparison with the post-LGM results, the grey lines in the right frames of Figs 3–5 are predictions of sea level change from the PGM to the LIG along each ridge. The broad trend in the prediction along the Mid-Atlantic Ridge across the earlier deglaciation is comparable to the post-LGM prediction (Fig. 2), however the predicted sea level variability along both the East Pacific Rise and Juan de Fuca Ridge is more muted. The former largely reflects the reduction in the gravitational migration of water away from North America during the penultimate deglaciation, which reduces the gradient in the far-field signal across the southeast Pacific. The latter reflects a reduction in both the gravitational signal and crustal subsidence at the periphery of the Cordilleran Ice Sheet, as well as the limited length of the ridge.

Fig. 7 shows predictions of sea level change from LGM to present at a series of additional ridges within the global ocean based on the 3-D GIA model. The change across the Galapagos Rift varies by less than ∼3 m across most of its extent, though it drops by ∼10 m as the ridge reaches Central America, where the prediction is influenced by ocean loading effects. Similar physics is active in the prediction for the Central Indian Ridge; the magnitude of the sea level change in this case drops significantly as the ridge enters the Gulf of Aden. The sea level change along the Chile Rise shows relatively muted variability except in the vicinity of the Patagonian Ice Sheet, where deformational and gravitational effects associated with ice unloading post-LGM are evident in a sharp reduction of ∼30 m. A reduction of greater magnitude is predicted on the western most section of the Gakkel Ridge, close to the Greenland Ice Sheet. Finally, the globe-encircling ridge system comprising the Southwest Indian Ridge, Southeast Indian Ridge, and Pacific-Antarctic Ridge is characterized by relatively minor variability, from ∼130–150 m.

3.3 Normal stress changes at mid-ocean ridges

As discussed in the Introduction, crustal stresses in submarine environments, including MORs, have commonly been associated—either implicitly or explicitly—with sea level (or meltwater) loading. The question as to whether the ice loading component of GIA drives significant crustal stresses in such environments has largely remained unanswered. To explore this issue, we compute normal stresses one computational grid element (14 km) below the surface and partitioned the results into contributions from ocean and ice loading. The calculation is technically challenging because the computed normal stresses are characterized by very large lateral gradients within the thin low-viscosity zone introduced in the Earth model as a proxy for mid ocean ridges and the calculation is sensitive to the details of the viscous structure adopted for this zone. Thus, raw profiles along MORs can exhibit large oscillations which we smooth using a simple moving window averaging. We note that the same approach was not necessary in the profiles in Figs 3–5 since the global sea level field is inherently smooth.

Fig. 8 shows profiles of (normalized) normal stress change from LGM to present day due to the ice and water loading components of the GIA process at the Mid-Atlantic Ridge, East Pacific Rise and the Juan de Fuca ridge. The results in the last case indicate—perhaps not surprisingly given the proximity of the ridge to the Cordilleran and Laurentide Ice Sheets—that crustal stress changes along the Juan de Fuca Ridge are strongly dominated by ice load induced deformation. The results for the Mid-Atlantic Ridge are more nuanced. In this case, the ice loading contribution only dominates a few 100 km from Iceland, whereas south of ∼50° latitude, the normal stress change associated with ocean loading is more significant. Some of the structure of the ice loading contribution below ∼60° N is clearly influenced by subsidence of the large-scale peripheral bulge of the Laurentian Ice Sheet (see Fig. 3). Finally, along the East Pacific Rise, the stresses associated with ocean loading dominate. The magnitudes are a factor of two different at the equator and an order of magnitude different at 30° S.

4. CONCLUSIONS

A long history of studies has demonstrated that sea level variations following the retreat of grounded ice sheets display significant geographic variability. These patterns reflect the complex gravitational, deformation and rotational effects of the GIA process. We have focused here on sea level changes across the last two deglaciation phases sampled by mid-ocean ridges that snake through the ocean basins, specifically the Mid-Atlantic Ridge, the Juan de Fuca Ridge and the East Pacific Rise. These ridges represent examples of slow, intermediate and fast spreading ridges, and they sample the geographic variability of post-glacial sea level in unique ways.

The predicted sea level change across the Mid-Atlantic Ridge shows the greatest variability, ranging from a maximum rise of ∼150 m in the South Atlantic to a sea level fall of ∼260 m in Iceland, and it is characterized by significant non-monotonicity as the ridge weaves in and out of the near-field gravitational and deformational (peripherical crustal subsidence) effects of GIA. The general trend in this prediction is similar across the penultimate and final deglaciation phases of the ice age. The predicted sea level trend across the East Pacific Rise is also significant, although the peak-to-peak range in the case of the penultimate glaciation (120–145 m) is half the range predicted across the final deglaciation event (100–150 m). Despite its proximity to the Cordilleran Ice Sheet, the variability in the predicted sea level change across the Juan de Fuca ridge is relatively limited. This is due to both a cancellation of GIA effects and the limited geographic extent of the ridge. Finally, the variability of sea level change at a series of other mid-ocean ridges (Fig. 7) is also relatively muted, though sharp gradients are evident in locations close to either continental margins (Galapagos Rift, Central Indian Ridge) and/or regions of ice mass flux (Gakkel Ridge, Chile Rise).

A link has been hypothesized, albeit with some contentious debate, between sea level fluctuations at mid-ocean ridges and various geodynamic processes including magma production and faulting, both of which have been connected to the development of ridge-perpendicular topographic variations, and hydrothermal fluxes. In addition, as noted in the Introduction, the first of these processes may provide a feedback in ice age climate. We have shown that any exploration of these links should also account for crustal stresses associated with ice loading which can extend well away from the location of ancient ice complexes. The results described herein demonstrate that the strength of these processes—if they have non-negligible impact on mid-ocean ridge dynamics—should exhibit the significant imprint of geographically variable glacial hydro-isostatic adjustment that has characterized the ice age world.

ACKNOWLEDGMENTS

This research was funded by Harvard University and the MacArthur Foundation. We are indebted to Giorgio Spada and an anonymous reviewer who provided numerous constructive comments on an earlier version of the manuscript. Prof. Spada suggested that we include ice load induced stresses on MORs.

AUTHOR CONTRIBUTIONS

Jackson Reilly (Conceptualization [equal], Formal analysis [equal], Investigation [equal], Methodology [equal], Visualization [equal], Writing - original draft [equal], Writing - review & editing [equal]), Konstantin Latychev (Investigation [supporting], Methodology [equal], Software [lead]), Sophie Coulson (Investigation [supporting], Methodology [supporting]), and Jerry X. Mitrovica (Conceptualization [equal], Formal analysis [equal], Funding acquisition [lead], Investigation [equal], Methodology [equal], Project administration [lead], Resources [lead], Software [lead], Supervision [lead], Validation [equal], Visualization [equal], Writing - original draft [equal], Writing - review & editing [equal]).

DATA AVAILABILITY

Results described in this paper are available on the public archive https://zenodo.org/records/14003374.

References