Imprint of relative sea level histories on Last Interglacial coral preservation

SUMMARY

Fossil corals are commonly used to reconstruct Last Interglacial (∼125 ka, LIG) sea level. Sea level reconstructions assume the water depth at which the coral lived, called the ‘relative water depth’. However, relative water depth varies in time and space due to coral reef growth in response to relative sea level (RSL) changes. RSL changes can also erode coral reefs, exposing older reef surfaces with different relative water depths. We use a simplified numerical model of coral evolution to investigate how sea level history systematically influences the preservation of corals in the Bahamas and western Australia, regions which house >100 LIG coral fossils. We construct global ice histories spanning the uncertainty of LIG global mean sea level (GMSL) and predict RSL with a glacial isostatic adjustment model. We then simulate coral evolution since 132 ka. We show that preserved elevations and relative water depths of modelled LIG corals are sensitive to the magnitude, timing and number of GMSL highstand(s). In our simulations, the influence of coral growth and erosion (i.e. the ‘growth effect’) can have an impact on RSL reconstructions that is comparable to glacial isostatic adjustment. Thus, without explicitly accounting for the growth effect, additional uncertainty is introduced into sea level reconstructions. Our results suggest the growth effect is most pronounced in western Australia due to Holocene erosion, but also plays a role in the Bahamas, where LIG RSL rose rapidly due to the collapsing peripheral bulge associated with Laurentide Ice Sheet retreat. Despite the coral model's simplicity, our study highlights the utility of process-based RSL reconstructions.

1 INTRODUCTION

Constraining the pattern and rate of global sea level change during the Last Interglacial (LIG; ∼128–116 ka; Rovere et al. 2016) offers insight into ice sheet sensitivity to temperature forcing, which underpins future sea level rise projections (e.g. DeConto et al. 2021). Fossil corals are widely used to reconstruct LIG sea level, as they grow close to the sea surface and can be accurately dated using U-series methods (Stirling & Andersen 2009; Chutcharavan & Dutton 2021). Sea level reconstructions based on fossil coral reefs estimate that global mean sea level (GMSL) during the LIG peaked at 5.5–9 m (Kopp et al. 2009; Dutton & Lambeck 2012) or less (Dyer et al. 2021) above the present day. However, significant debate surrounds the timing, duration, number and magnitude of the highstand(s) (e.g. Stirling et al. 1995, 1998; Neumann & Hearty 1996; Thompson et al. 2011; Kopp et al. 2013; Barlow et al. 2018; Dyer et al. 2021).

Coral growth is highly dependent on the coral's depth below the sea surface, and thus is influenced by rates of sea level change. Although growth rates vary by species, coral growth is dependent on light availability (Hopley 2011) such that the rate of vertical coral accretion can, for many species, be approximated as a decreasing logarithmic function of depth (Kleypas 1997). The maximum coral reef accretion rate is species dependent, but often occurs around 5 m depth, where wave stress is low and light availability is high (Woodroffe & Webster 2014). As sea level rises, a coral may continue growing to keep up with sea level (Hopley 2011), but will die if the long-term rate of sea level rise exceeds the long-term vertical accretion rate. When sea level falls, a coral may be eroded due to wave action (Toomey et al. 2013).

Sea level reconstructions using coral reefs rely on the relative water depth (also referred to as the indicative meaning), which is an assumption of how far below the mean lower low water the coral reef grew. This value is often based on modern coral depth distributions (e.g. Hibbert et al. 2016). Typically, the mean relative water depth and a corresponding uncertainty is added to the present-day coral fossil elevation to reconstruct relative sea level at the time of interest, when the coral reef was alive (e.g. Rovere et al. 2016). However, this approach assumes a relative water depth for a given reef that is constant through time and space. It thus fails to account for how coral growth may influence the relative water depth in response to spatially varying sea level change, and how coral erosion may expose coral reef surfaces with older ages and different relative water depths. Prior studies have used LIG coral erosional surfaces and reef growth patterns to assess sea level highstands or trends of sea level rise and fall at a given location (e.g. O'Leary et al. 2013; Skrivanek et al. 2018; de Gelder et al. 2022). Forward numerical models of coral reefs have also been used to better understand the impact of sea level oscillations on reef growth and preservation at one or two individual locations (Pastier et al. 2019; de Gelder et al. 2022; Boyden et al. 2023). Nevertheless, the integrated coral growth and erosion history over an entire glacial cycle and across locations has not yet been systematically considered.

Because corals are sensitive to sea level change, understanding coral growth, erosion, and ultimate preservation resulting from a range of possible LIG global ice volume scenarios can provide insight into a coral's ability to record local LIG sea level. In this study, we simulate coral growth and erosion in response to sea level change using a simple coral parametrization to model the formation and preservation of LIG corals in the Bahamas and western Australia, regions which house a large number (n > 100) of LIG fossil corals and experienced different relative sea level (RSL) histories over the last glacial cycle given their respective location in relation to past former ice sheets (Dendy et al. 2017; Fig. 1). Since the history of global sea level change is uncertain, we explore a broad range of possible scenarios over the LIG by varying the magnitude, timing, and pattern of peak GMSL. Our goal is not to robustly constrain LIG sea level or interpret specific fossil coral reef records; rather, we seek to discover how coral growth and erosion interplays with relative sea level change across the Bahamas and western Australia to produce present-day preserved LIG coral fossil elevations. Quantifying this relationship may provide additional insight into the ability of preserved present-day corals to capture local LIG sea level, as different patterns of local LIG sea level change may influence a coral's depth below the sea surface, as well as the uncertainty associated with coral-based sea level reconstructions.

2 METHODS

2.1 Modelling coral growth and erosion in response to sea level change

Our goal is to isolate the effects of sea level change on coral growth patterns. Thus, we choose a simplified coral model with growth and erosion rates based on averaged accretion rates during the LIG of corals in the Pacific, Atlantic, and Indian Oceans (Toomey et al. 2013). The corals in the LIG data set (Fig. 1a, yellow dots) experienced a complex interplay of processes besides sea level change, including interspecies competition, substrate availability, observation bias, spatially heterogenous wave regimes and detailed bathymetry (Woodroffe & Webster 2014), which are not accounted for in our generalized coral model.

The model we use is based on the mathematical relationship between coral surface elevation and vertical accretion rate, first suggested by Bosscher & Schlager (1992) and updated by Toomey et al. (2013) to include the effects of tidal depth. For each time-step |$(\Delta t$| = 10 yr), the model calculates the vertical accretion rate of existing live corals in response to local sea level change |$\frac{{{\rm d}S}}{{{\rm d}t}}$|⁠, which is comprised of a GMSL change signal and a GIA signal. The vertical growth or erosion rate (⁠|$\frac{{{\rm d}Z}}{{{\rm d}t}}$|⁠) is described by:

Here Z is the coral elevation [m|$];\ {G}_{\rm max}$| is the maximum vertical accretion rate [m yr⁻¹]; the constant, 4.4, is an averaged |$\frac{{{I}_0}}{{{I}_k}}$| where |${I}_0$| is surface light intensity and |${I}_k$| is saturating light intensity [unitless]; D is tidal depth (mean tidal range) [m]; W is erosion [m yr⁻¹] and |${s}_ + $| and |${s}_ - $| are the positive (above the sea surface) and negative (below the sea surface) boundaries for the surf zone. The values for the parameters adopted in this study are shown in Table 1.

We initialize coral models at 132 ka (∼2–7 ka before the peak LIG GMSL). This assumes that sea surface temperatures in western Australia and the Bahamas at this time were similar to present, allowing for coral growth (Hoffman et al. 2017). We simulate a single reef in each grid cell with observed LIG corals (0.35° |$\times $| 0.35° resolution). Simulated corals are placed at the tidal depth associated with sea level at the initiation time (see Table 1). Coral growth occurs when the coral is underwater (Z < 0), and corals stop growing when their elevation falls below a threshold called the limiting depth (Z < Z_L) or when the rate of vertical accretion is less than the rate of sea level rise (⁠|$\frac{{{\rm d}Z}}{{{\rm d}t}}$| < |$\frac{{{\rm d}S}}{{{\rm d}t}}$|⁠; Kleypas 1997; Toomey et al. 2013). If corals are exposed subaerially (Z > 0), corals stop growing, and coral growth is re-initiated when corals are re-submerged. After each simulation, the age of the modelled corals are calculated as the time at which the exposed coral surface was last in growth position (i.e. the most recent time when Z < 0). The RSL recorded by the coral (i.e. RSL at the time the coral was most recently growing) is then calculated as the RSL at that time (see Fig. 2).

We then identify the impacts of dynamic coral growth and GIA on modelled present-day elevations of LIG coral fossils. We define the growth effect as the RSL recorded by the coral subtracted from the present-day coral elevation, which is equivalent to the water depth at which the coral grew. We define the GIA effect as the GMSL subtracted from the RSL recorded by the coral. The growth effect thus includes coral growth and erosion in response to sea level changes, and the GIA effect does not. The growth effect is always negative because our model does not permit coral growth above the sea surface. Thus, the present-day elevation of the coral is the RSL recorded by the coral plus the (negative) growth effect. The growth effect has a maximum value of tidal depth (−1.15 m in the Bahamas, −1.50 m in western Australia); a growth effect at this value indicates that the coral was fully caught up to LIG RSL at the time of its death. A more negative growth effect indicates that the coral was not caught up with RSL at the time of its death, and tidal depth cannot be used to reconstruct RSL. Fig. 2 shows a schematic of coral growth and erosion through time, and resulting growth and GIA effects.

2.2 Glacial isostatic adjustment model and global mean sea level over the last glacial cycle

To predict the relative sea level history for each LIG coral site in the Bahamas and western Australia, we first generate a set of 900 GMSL histories over the last glacial cycle. This allows us to explore impacts of GMSL peak type, peak timing, and peak magnitude on the relative sizes of the growth and GIA effects in the Bahamas and western Australia. We assume that the penultimate glacial maximum had the same GMSL as during the Last Glacial Maximum (−130 m), and that the penultimate deglaciation lasted ∼20 kyr (putting the penultimate glacial maximum at 150 ka), although there is evidence that the penultimate deglaciation was shorter than 13 kyr (e.g. P. U. Clark et al. 2020). After pinning the penultimate glacial maximum at 150 ka with a GMSL of −130 m, we sample a range of possible GMSL histories during the LIG, based on previously suggested patterns and magnitudes of GMSL change (e.g. Neumann & Macintyre 1985; Stirling et al. 1998; Hearty et al. 2007; Rohling et al. 2008; Blanchon et al. 2009; Thompson et al. 2011; Kopp et al. 2013; Dechnik et al. 2017; Skrivanek et al. 2018). The peak GMSL magnitude is randomly sampled between 0 and 10 m, where the latter is the 15 per cent probability exceedance value of the LIG GMSL highstand estimated by Kopp et al. (2009). and the lowstands are sampled between −2 m and the lowest sampled GMSL peak for a given run.

We also vary the pattern of GMSL during the LIG according to three scenarios: a single GMSL highstand referred to as ‘single peak’ (Fig. 3b), a delayed ascending GMSL highstand preceded by a 2 kyr period of stability called ‘ascending peak’ (Fig. 3c), and a double GMSL highstand separated by 2 kyr of lower GMSL, called ‘oscillating peak’ (Fig. 3d). In addition, we vary the timing of the GMSL highstand, which modulates the rate of sea level change preceding and following the LIG GMSL highstand (Table 2). Across the LIG (between 130 and 115 ka), the rate of GMSL change varies from −7.7 to 16.8 m kyr⁻¹.

We sample a range of estimates during MIS 5a–d to populate GMSL history through 70 ka (Fig. 3a; Cutler et al. 2003; Creveling et al. 2017). Sea level falls at a rapid, intermediate, or slow rate to −60, −40 or −20 m, respectively, during MIS 5d (110 ka). To account for potential erosion or re-submergence of LIG corals, we vary the magnitude of subsequent GMSL highstands during MIS 5c (100 ka), although our subsequent analysis only considers modelled reefs aged 132–115 ka. The GMSL highstand during MIS 5c is either high-, intermediate- or low-magnitude, with a value of 0, −10 or −20 m, respectively. The sampled GMSL values for MIS 5d and MIS 5c are repeated to fill in the loading history during MIS 5b (90 ka) and MIS 5a (80 ka), respectively (Fig. 3a; dashed black lines). We use shape-preserving piecewise cubic interpolation to fill the GMSL curve between sampled points from 150 to 70 ka. From 70 ka until modern day, we adopt the GMSL history associated with the ICE-5 G model (Peltier 2004). While the corals are exposed for much of this period, including this period allows us to account for potential erosion or re-submergence of corals during the late Holocene (∼8 ka to present). This is especially important for sites in western Australia, which experienced a local highstand at the end of the last deglaciation phase rather than at present-day. Although previous studies have shown that the GMSL history associated with ICE-5 G underestimates sea level during MIS 3 (Pico et al. 2016), our modelled coral distributions are insensitive to GMSL oscillations below ∼−20 m, the limiting depth in the coral model (see Table 1).

We then produce a set of ice histories corresponding to the ensemble of 900 GMSL histories (Fig. 3). The ice histories are based on ICE-5 G (Peltier 2004), and assume that ice geometry prior to 26 ka is identical to ice geometry post-26 ka for the same GMSL value. During the LIG, when GMSL values exceed present-day GMSL, we uniformly reduce the thickness of the Antarctic and Greenland ice sheet to produce excess melt consistent with the sampled GMSL values. There are some limitations to this approach; for example, there is evidence that ice volumes over North America and Eurasia differed between the last glacial maximum and penultimate glacial maximum (Colleoni et al. 2016; Rohling et al. 2017), and subsequent research has also refined the deglacial ice thickness history over North America, Northwestern Eurasia and Antarctica, resulting in updated global ice histories (e.g. Roy & Peltier 2015, 2018). However, the focus of our analysis is sea level prior to 26 ka, a period characterized by large uncertainties in ice extent and volume (e.g. Dalton et al. 2022), and therefore we do not focus on the impact of different ice geometries on reef development.

We reconstruct RSL at each site in the Bahamas and western Australia by performing glacial isostatic adjustment (GIA) simulations with each of the 900 ice loading histories. Our calculations are based on the theory and pseudo-spectral algorithm described by Kendall et al. (2005) with a spherical harmonic truncation at degree 256. These calculations adopt a Maxwell rheology that is incompressible in the fluid limit and include both the impact of load-induced Earth rotation changes on sea level (Milne & Mitrovica 1996, 1998) and evolving shorelines, where the latter incorporates the evolution of grounded, marine-based ice (e.g. Milne et al. 1999; Kendall et al. 2005).

In addition to the history of global ice cover, our predictions require regional models for Earth's viscoelastic structure. The Earth model for the Bahamas is characterized by a lithospheric thickness of 96 km, an upper mantle viscosity of |$5 \times {10}^{20}$| Pa s and a lower mantle viscosity of 4 |$\times {10}^{21}$| Pa s, which is consistent with MIS 5a and 5e sea level data in the Western North Atlantic and other estimates for Earth structure in this region (Potter & Lambeck 2004; Creveling et al. 2017), although higher viscosity estimates in the upper and lower mantle have been obtained using Holocene data (e.g. Milne & Peros 2013). The earth model for western Australia is characterized by a lithospheric thickness of 71 km, an upper mantle viscosity of |$5 \times {10}^{20}$| Pa s and a lower mantle viscosity of |${10}^{22}$| Pa s, which is consistent with late Holocene sea level records in Australia (O'Leary et al. 2013).

To focus on the influence of LIG GMSL histories on coral preservation, we limit our study to a single Earth model for the Bahamas and a single Earth model for western Australia and a single ice history for the penultimate deglaciation extending to 150 ka. Although a suite of 900 sea level histories are produced, we present results only from histories that generate modeled LIG reef fossil surfaces (see Section 3.4).

2.3 Last Interglacial fossil corals in the Bahamas and western Australia

Fig. 1(a) shows a global database of fossil corals (Hibbert et al. 2016). From this database, we compile LIG fossil coral elevations from the Bahamas (77.3–73.7°W and 21.0–26.4°N) and western Australia (113.1–115.5°E and 34.2–21.8°S), two regions characterized by a large sample size (n > 100) of LIG fossil corals that have U-series ages between 115 and 130 ka. Although these regions are widely recognized as passive margin settings, tectonic activity has been observed in western Australia (e.g. Sandstrom et al. 2020). Corals are excluded from the database if they are explicitly listed as not in growth position in the original reference. We note that many of the coral samples do not pass open-system screening requirements (Hibbert et al. 2016), which can bias the ages. In the Bahamas, 155 samples are located among three sites, and in western Australia 102 samples are located among 15 sites (Figs 1b and c). To reduce spatial bias in the observed distributions resulting from repeated sampling from the same site, we adopt a sampling frequency based on kriging weights. We use ordinary kriging and a theoretical exponential variogram to calculate these weights, and obtain a sampling frequency by setting the standard deviation of the weights to 20 per cent of the mean, sampling approximately 10 000 times in each region. Coral elevation distributions across the Bahamas and western Australia after sampling are shown in Figs 4(a) and (b), respectively. These distributions are markedly different; for instance, the distribution in western Australia is characterized by a wider range and higher mean value than that in the Bahamas.

3 RESULTS AND DISCUSSION

3.1 Relative sea level patterns in the Bahamas and western Australia

The Bahamas and western Australia experienced different RSL histories during the LIG due to GIA (Fig. 1). If global mean sea level was stable across the LIG, Western Australia would experience an early RSL highstand followed by an RSL fall due to the far-field effects of continental levering and ocean syphoning (Dutton & Lambeck 2012; O'Leary et al. 2013; Fig. 5). In contrast, RSL would rise throughout the LIG in the Bahamas due to the ongoing collapse of the peripheral bulge of the Laurentide Ice Sheet (Potter & Lambeck 2004; Fig. 5). Proximity to the Laurentide Ice Sheet also governs the spatial pattern of RSL within the Bahamas; the north-south gradient in Fig. 1(b) reflects the difference in the level of isostatic equilibrium at 125 ka across the bulge relative to the present day (Potter & Lambeck 2004; Creveling et al. 2017; Dyer et al. 2021). In contrast, RSL in western Australia (Fig. 1c) does not show a significant gradient between sites because GIA across the coastline reflects an ocean loading signal that is roughly shoreline-parallel (Nakada & Lambeck 1989; O'Leary et al. 2013). Differences between coral elevation distributions in the two regions (Fig. 4) may in part be explained by differences in their RSL history (Fig. 5).

3.2 Modelling preserved coral elevations at a single site

We begin with a case study to explore coral growth and preservation at two sites, Abaco Island in the Bahamas (star, Fig. 1b) and Houtman Abrolhos Islands in western Australia (star, Fig. 1c). At each location, we force our coral model with RSL associated with two different GMSL curves. The GMSL histories have the same highstand (∼9 m) and peak type (ascending peak), but different timing (early, Figs 6a, c, e; and intermediate, Figs 6b, d, f) and different magnitudes of initial peaks (higher first peak in Figs 6a, c, e; lower initial peak in Figs 6b, d, f). Hereafter, the GMSL histories are referred to as ‘early peak timing GMSL’ and ‘intermediate peak timing GMSL’. For each simulation, our model produces a present-day coral elevation (Figs 6a–d, red dot), which can be corrected to reflect LIG sea level (Figs 6a–d, blue dot). The modelled elevation at present day is the elevation one would expect the top of the LIG reef to be observed today if it had experienced the modelled history of growth and erosion. The difference in coral elevation in Fig. 6(b) between the top of the reef at 112 ka and the present-day reflects subsequent erosion of the reef during MIS 5a–d.

We compare modelled coral elevations using two RSL histories generated by the early peak timing GMSL (Fig. 6a) and intermediate peak timing GMSL (Fig. 6d). In our model, coral growth occurs when corals are below the sea surface and the rate of coral accretion exceeds the rate of sea level change (e.g. Figs 6b, c, e, f). Coral erosion occurs when (1) sea level falls (e.g. Figs 6b, c, e, f), (2) corals are in the surf zone during a period of rapid sea level rise that outpaces the maximum rate of coral accretion and corals cannot keep up (e.g. Figs 6b, e, f) or (3) when post-LIG sea level is high enough to erode LIG corals (e.g. Fig. 6e). Coral elevation is static when (1) the rate of sea level rise exceeds the coral growth rate and thus corals cannot grow, but corals are not in the surf zone and thus cannot be eroded (e.g. Figs 6b, e, f), (2) corals have caught up with sea level and can no longer grow towards the sea surface or (3) corals are subaerially exposed, but are not in the surf zone and thus cannot be eroded (e.g. Figs 6b, c, e, f). On Fig. 6, periods of coral erosion are shaded in grey, periods of coral growth are shaded in yellow, and periods of static coral elevation are unshaded.

The results show that RSL histories with similar GMSL highstands can produce modelled coral fossils with >1 m difference in present-day elevation due to GIA and the dynamic coral response to sea level change (growth effect). In the Bahamas (Figs 6b and e), the two GMSL histories produce RSL peaks that differ by 0.6 m (10.3 m in Fig. 6b; 10.9 m in Fig. 6e). However, the preserved coral fossils both record RSL of ∼10 m (9.9 m in Fig. 6b, 10.2 m in Fig. 6e). Thus, the RSL peak is better captured by the coral experiencing the early ascending peak GMSL (Fig. 6b) compared with the coral experiencing the intermediate ascending peak GMSL (Fig. 6e). This effect is largely due to differences in post-LIG RSL; in the intermediate ascending peak GMSL, subsequent erosion of LIG corals produces a large growth effect (Fig. 6e). We also note that while the present-day coral elevations differ by ∼1 m (7.0 m in Fig. 6b; 6.0 m in Fig. 6e), both corals record the same magnitude of RSL, suggesting that different present-day coral elevations are not necessarily indicative of different recorded sea levels.

In western Australia, at Houtman Abrolhos Islands (Figs 6c and f), the RSL peaks are equivalent (7.5 m) across the two GMSL histories. This RSL peak is well-captured by the coral experiencing the intermediate ascending peak GMSL (Fig. 6f; RSL recorded by the coral is 7.0 m), but is not well-captured by the coral experiencing the early ascending peak GMSL (Fig. 6c; RSL recorded by the coral is 5.7 m). With the early ascending peak GMSL, RSL falls more gradually after the peak GMSL, leading to a longer period of coral erosion and a lower present-day coral elevation. This indicates that even when no post-LIG erosion occurs, the same RSL peak can produce different present-day coral elevations depending on rates of LIG sea level change. Thus, the growth effect has the potential to introduce additional uncertainty into sea level reconstructions since the corrections to reconstruct LIG sea level using our process-based coral model depend on the GMSL histories.

3.3 Modelling preserved coral elevations across sites

The magnitude of the growth effect (i.e. the difference between blue and red dot in Figs 6b, c, e and f) varies with location, and differs between the near-field (Bahamas) and the far-field (western Australia). We next examine simulated coral elevations across the Bahamas and western Australia produced by the same two GMSL histories. The early ascending peak GMSL (Fig. 6a) and intermediate ascending peak GMSL (Fig. 6b) produce distinct predicted LIG coral elevations, GIA effects, and growth effects across all sites in the Bahamas (Figs 7a–c) and western Australia (Figs 7d–f).

In the Bahamas, coral elevations in the Bahamas simulated by intermediate ascending peak GMSL are lower than those simulated by the early ascending peak GMSL because the intermediate ascending peak GMSL has higher sea level during MIS 5a–d, leading to post-LIG erosion of some LIG corals (corals that eroded after LIG are indicated by triangles on Fig. 7a). Generally, corals that eroded after the LIG have lower elevations, although this is not always the case (e.g. green squares in Fig. 7a). The growth effect varies across sites (e.g. Fig. 7c), most notably due to erosion; corals that experienced post-LIG erosion are associated with significantly larger-magnitude growth effects. This demonstrates that gradients of coral elevations across sites—which have been assumed to be representative of GIA gradients (e.g. Potter & Lambeck 2004)—can be produced by an interplay of GIA and coral growth and erosion, although we note that on large spatial scales the GIA effect will dominate the growth effect. This is particularly important to consider when coral elevations are interpreted without consideration (or elimination) of corals that show erosion.

In western Australia, the growth effect varies more dramatically within a single simulation (∼4.5 m range in growth effect when forced with the intermediate ascending peak GMSL). As in the Bahamas, the corals that experienced post-LIG erosion are associated with significantly larger-magnitude growth effects (Fig. 7f). We also note that the growth effect across the two GMSL histories varies more widely in western Australia due to more significant post-LIG erosion at some sites when the corals are forced with the intermediate ascending peak GMSL.

3.4 Modelled coral elevation distributions using a range of LIG global ice volume histories

Next, we force our coral model at each LIG coral site with the relative sea level histories associated with the 900 GMSL histories (Fig. 2). To extract LIG corals from our simulation, we filter for reefs with exposed surfaces last growing between 132–115 ka. Of 900 simulations, 295 produce modelled LIG reefs at all sites in the Bahamas and 11 of the 15 sites in western Australia where observed LIG corals exist. The ranges of modelled coral elevations produced across all these runs at each site in the Bahamas and across western Australia are shown in Fig. 8(a) (solid bars). Ranges of observed coral elevations are shown for comparison (Fig. 8a, dotted bars).

We examine how different GMSL histories produce different modelled coral elevations. We present subsets of the full 295 simulations in Figs 8(b)–(d) to isolate the impacts of LIG peak timing, peak magnitude, and peak type on the modelled coral elevations. To begin, we consider modelled corals produced by a range of highstand values, controlling for type of GMSL peak and timing (ascending GMSL peaks and early GMSL peak timing; Fig. 8b). As the highstand value increases, the upper bounds on the coral elevation ranges increase at all sites. We next consider how modelled coral elevation varies with peak type (Fig. 8c) and timing (Fig. 8d) by holding the GMSL highstand constant at ∼7.8 m. The range of modelled coral elevations in western Australia is similar across these subsets (Fig. 8; blue). However, the range of modelled coral elevations in the Bahamas is highly sensitive to variations in peak type and peak timing; single peaks and early peaks produce a lower range of modelled coral elevations at all three Bahamian sites than ascending or oscillating GMSL peaks or GMSL peaks with intermediate or late timing (Figs 8b–d).

3.5 Comparing the impact of glacial isostatic adjustment (GIA) and coral dynamics

Different model runs produce different magnitudes of the growth effect (growth effect varies between ∼–8.5 m to tidal depth; Figs 8h–j). A wider range of the growth effect implies more uncertainty in the depth at which corals at that site grew below the sea surface. Because the growth effect is not constant across runs, even after controlling for peak type, peak magnitude and/or peak timing, RSL reconstructions that rely on coral fossils may be subject to additional uncertainty resulting from an unknown growth effect. For example, if a growth effect of tidal depth is assumed to be constant across sites (implying that corals were fully caught up to RSL when they died), reconstructions could underestimate RSL by up to ∼7 m. The growth effect (and thus the potential to underestimate RSL) tends to be higher in magnitude for corals that experienced post-LIG erosion (triangles, Fig. 8h–j); however, even for corals that did not experience post-LIG erosion (squares, Figs 8h–j), reconstructions could underestimate RSL by up to ∼4.5 m.

In our model, higher-magnitude GMSL peaks increase the upper bound on coral elevations; as GMSL highstands increase, corals respond by growing to higher elevations. More prolonged peaks (ascending or oscillating GMSL peaks) increase the upper and lower limits of coral elevation at every site in the Bahamas, allowing more time for the corals to catch up to sea level than the shorter single peaks. The modelled elevation range at the Bahamas is sensitive to relative sea level gradients along the peripheral bulge of the North American ice sheet (Fig. 1b), illustrated by the differences in the GIA effect across Abaco Island, San Salvador Island, and Great Inagua Island (Fig 8e–g). GIA has a ∼–5 to +5 m effect on modelled coral elevations in the Bahamas (Figs 8e–g). The growth effect is mostly ∼–3 m to tidal depth (–1.15 m) in the Bahamas, with a wider spread of growth effect values at the northernmost site (Abaco Island, down to −5.5 m) due to its location on the peripheral bulge. Abaco Island is characterized by the highest GIA-induced sea level rise, and the reefs struggle to keep up with the more rapid sea level rise.

Because LIG RSL in western Australia is characterized by an early peak followed by a relative sea level fall (Fig. 5), corals are more likely to keep up with sea level with a slower deglaciation (and later GMSL peak timing). The GIA effect here varies between ∼−3 to +3 m. For each GMSL subset, the range of the GIA effect across western Australia is wider than the range of the GIA effect across the Bahamas, despite the larger GIA gradient across sites in the Bahamas (Figs 1b, cf. Fig. 1c). This feature is due to post-LIG erosion in western Australia. There is wider age distribution in western Australia because Holocene RSL rose above present-day sea level (Fig. 5), which erodes the most-recent coral, and creates LIG coral surfaces with older ages than in the Bahamas (Fig. 9). The corals in western Australia sample RSL across a greater duration of the LIG, and are thus characterized by a wider range of GIA corrections.

This simulated Holocene erosion also results in the wider range of the growth effect magnitude in western Australia, since corals can be eroded to older ages when the reef had not yet caught up to sea level. These results indicate that the growth effect varies with the cumulative history of GIA effects, and is sensitive to both RSL during the LIG and subsequent highstands which may have eroded the preserved corals.

While the upper bound on modelled coral elevation range in western Australia is substantially lower than the observed coral elevation range (11.1 m for the observed range and 6.7 for the modelled range), only 9 out of 102 samples in the data set of observed elevations are >6.6 m, and these fall more than two standard deviations away from the mean. Five of the nine samples are located at Cape Cuvier and Cape Range, where several prior studies have argued there to be Quaternary tectonic activity (e.g. Condon et al. 1955; Clark 2010; Whitney & Hengesh 2015; Sandstrom et al. 2020). If we remove elevation data in western Australia that fall more than 2 standard deviations away from the mean, the upper bound on the range of observed coral fossil elevations is 6.6 m, similar to the modelled range (black triangle; Fig. 8a). There are also two locations (San Salvador Island and Vlaming Head in western Australia) where the lower bound on observed coral elevations does not overlap with the range of modelled coral elevations. This may be due to our assumption that corals are caught up with sea level at the time of their initialization.

None of the simulations produce LIG reefs at four of the sites in western Australia (Burney Point, 2 sites at Shark Bay and Leander Point), where 19 coral fossil data points are observed. This is likely due to erosion of coral fossils in our model during the higher Holocene RSL in western Australia. Erosion is a testable feature in our model; for example, past studies have suggested that Shark Bay coral communities have high potential for ongoing erosion (O'Leary et al. 2008). Because of the simplifying assumptions we make in constructing our coral model, our modelled coral results should not be quantitatively compared with observations. However, locations of post-LIG coral erosion in our model can be tuned to locations of eroded corals in the observational record, which would allow for a qualitative comparison between our coral model results and field observations.

3.6 Modelling limitations

To isolate the impacts of GMSL on coral growth and erosion, our study makes a set of simplifying assumptions. First, our coral model uses constant maximum accretion and erosion rates, although these values are unlikely to be the same across coral species. Our model also does not include lateral accretion or allow for different wave regimes depending on local topography, which would influence reef growth and preservation (Woodroffe & Webster 2014). An example of a more complex reef evolution model which includes topographic and sediment effects can be found in Salles et al. (2018). Including these effects would likely spread the elevation distribution of modelled preserved coral fossils, and potentially would allow for direct comparisons between modelled coral fossil elevations and observed coral fossil elevations.

Secondly, our GIA simulations adopt a single ice distribution configuration and a single Earth model for each region. Varying the region and timing of excess ice melt during the LIG, including generating ice histories with different ice geometries, would impact the shape of modelled coral distributions given the sensitivity of our coral model to RSL changes over the LIG. For example, a collapse of the West Antarctic Ice Sheet at different times during the LIG would influence RSL histories in western Australia and the Bahamas, and these would produce distinctive modelled coral distributions (Hay et al. 2014). In addition to being sensitive to the melt source during the LIG, these results are also sensitive to the ice distribution prior to the LIG. RSL at the Bahamas, and therefore the distribution of coral elevation at the Bahamas, is highly sensitive to the size of the MIS 6 Laurentide ice sheet and its speed of deglaciation (Dendy et al. 2017). Accounting for different temporal and spatial patterns of ice sheet melt across the LIG on coral elevation distributions, in addition to different viscoelastic Earth structures, should be the subject of future study.

5 CONCLUSIONS

By modelling coral growth and erosion in response to sea level change, we demonstrate that the elevations of preserved LIG corals in the near- and far-field are sensitive to the magnitude, rate and timing of global mean sea level change over the LIG. In particular, we perform a suite of GIA simulations spanning a wide range of possible LIG sea level histories and simulate the evolution of LIG corals in response to the relative sea level change in the Bahamas and western Australia. This allows us to systematically isolate the influence of GIA on modelled coral elevations, which is distinct from the influence of dynamic coral growth and erosion.

Although this work makes a series of simplifying assumptions, our findings illustrate that coral growth and erosion in response to local sea level change influences the coral's preserved age, present-day elevation, and depth below the sea surface at its time of most recent growth (which we call the growth effect). Across our simulations, the growth effect has a maximum magnitude of 5.4 m in corals that did not experience post-LIG erosion and can be even larger (8.9 m) in corals that did experience post-LIG erosion, although the mean growth effect is significantly smaller (approximately 1.6 m in the Bahamas and 2 m in western Australia, for both eroded and non-eroded corals). The dynamic response of coral growth to sea level change introduces uncertainty into RSL reconstructions, with the magnitude of this uncertainty dependent on chosen parameter values, location and sea level history.

To constrain the true magnitude of the growth effect, future work could invoke an inverse coral model to reconstruct RSL using the present-day distribution of LIG corals. Field observations of post-LIG coral erosion could also inform estimates of the growth effect. In addition, an improved coral model which incorporates more complex processes, including species-specific growth and erosion patterns and lateral accretion, may be used to explore a range of GMSL scenarios and melt patterns, producing modelled elevation distributions which can be compared to observed elevation distributions in both the Bahamas and western Australia. Such exercises could demonstrate the potential of coral modelling as an innovative tool to reduce uncertainty in coral-based LIG RSL reconstructions, given that only certain GMSL histories will produce good fits between the observed and modelled elevations. One advantage to such an approach is that comparing modelled and observed coral elevation distributions does not require good age control beyond an association with the LIG. While corals tend to have good age constraints, they can be plagued by open-system behaviour. Furthermore, other LIG sea level indicators such as marine terraces and erosional notches have significantly larger age uncertainties and applying a similar process-based approach may also provide new insight from these data sets.

ACKNOWLEDGEMENTS

JA acknowledges support from NSF award MGG 18–41888. Funding was provided by Harvard University andNSF grant OCE-1702684 (JXM). TP acknowledges support from Harvard University, NSF EAR Postdoctoral Fellowship and UC President's Postdoc Program Fellowship. RCS, TP, JXM, JA and PJH designed the study and wrote the paper. RCS performed the modelling and analysis. We thank the editor, Yucheng Lin and an anonymous reviewer for their comments and improvements to the paper.

DATA AVAILABILITY

The Last Intergacial coral fossil elevation database used in this study (see Figs 1 and 4) is archived by British Oceanographic Data Centre (BODC; www.bodc.ac.uk) doi: 10.5285/32056c4c-fef8-29c4-e053-6c86abc06cd4 (Hibbert et al. 2016). Global mean sea level histories used in this study and a MATLAB version of the coral model will be made available at https://github.com/becca-cs/lig_coral.

References