Fossil Java Sea corals record Laurentide ice sheet disappearance

Outside regions of major tectonics, sea-level change during glacial cycles is driven by the mass exchange between grounded ice sheets and the oceans, and its associated planetary deformation. The isostatic signature of this mass exchange varies geographically, according to whether the ice or water load is dominant, and locally, according to rheological parameters that are depth-dependent and laterally variable (Mitrovica and Peltier, 1991; Lambeck, 1993; Lambeck and Chappell, 2001; Mitrovica and Milne, 2002; Woodroffe et al., 2012). Predictive models of past and ongoing isostatic responses are important to constrain key processes of future relative sea-level (RSL) rise due to global warming. Such models rely on information about the ice-sheet evolution and the planet’s rheology.

Glaciological chronostratigraphic constraints provide information on the evolution of large ice masses such as the vast North American Laurentide ice sheet (LIS). Typically based on ¹⁴C dating of organic material beneath glacial till deposits, LIS deglaciation was identified to have occurred between 8 and 6 ka (Dyke and Prest, 1987; Dyke et al., 2002). Glaciological constraints have been adopted, for example, by Peltier (2004), Peltier et al. (2015), and Lambeck et al. (2017) to reconstruct the deglaciation chronology in their isostatic models. Due to the unknown time that elapsed between deposition of the dated material and ice retreat, however, ¹⁴C data provide minimum-limiting information only. Later studies circumvented this problem through cosmogenic ¹⁰Be dates and constrained LIS melting to 7.1–6.3 ka (Carlson et al., 2008; Ullman et al., 2016).

Geological observations of past RSL change provide further information on ice dynamics to calibrate isostatic models. Importantly, the inversion of such evidence from regions distant to the former ice sheets (i.e., far field) provides data to quantify both the parameters of the ice-sheet evolution that are missing from the glaciological database and the rheological parameters that define the response to ice-water surface loading (Lambeck, 1993). In the global model of glacial isostatic adjustment applied in this work, inversion of far-field RSL data, and its associated uncertainties, was interpreted to indicate the end of significant Northern Hemisphere deglaciation at 6.8 ka (Lambeck et al., 2014).

This study presents a synthesis of RSL data inferred from mid-Holocene coral microatolls at the tectonically stable north coast of Java, Indonesia. Upward growth of nonponded microatolls is constrained by low water levels, rendering their upper surfaces as precise indicators for RSL change (Scoffin and Stoddart, 1978; Smithers and Woodroffe, 2000). Our results provide new constraints with which to calibrate predictive models of isostatic response. Furthermore, results are of societal importance because instrumental and proxy records of RSL change are short in Indonesia, yet the impacts of future RSL rise on communities in this region will likely be profound.

Pulau Panjang is a sand cay on top of a coral reef platform in the marginal Java Sea 2 km offshore Jepara at the north coast of central Java, Indonesia (Fig. 1). The windward northwestern margin of the intertidal reef flat has a rubble sheet, and the vegetated island is located on the central platform. Living coral occurs along the windward reef margin. The tidal cycle is mixed semi-diurnal with a maximum tidal range of 1.4 m.

The Java Sea is a well-suited far-field region for RSL reconstructions: First, it is located outside the storm belt. Consequently, large corals interpreted as RSL indicators in this study have likely not been relocated by storm waves. Second, the Java Sea forms the southeastern margin of the tectonically stable Sunda Shelf. Jepara is located more than 400 km north of the Sunda megathrust, and so tectonic uplift is unlikely (Sobolev et al., 2007). Around the study site, only a few deep-seated earthquakes have been recorded, all of minor magnitude and without vertical displacement effects at the surface (Simons et al., 2007; see also Supplemental Material¹). Third, only one previous RSL reconstruction exists nearby, at Teluk Awur (Azmy et al., 2010), and after data standardization, this record provides evidence of higher RSL (+~1.3 m) that remained constant (within the observational accuracies) between 6.7 and 5.3 ka (Mann et al., 2019).

Our RSL reconstruction followed a standardized approach with verifiable age and elevation uncertainties (Shennan et al., 2015; Khan et al., 2019). We constrained mean sea level (MSL) at Pulau Panjang to 0.76 ± 0.12 m (2σ) above the World Geodetic System 1984 (WGS84) ellipsoid based on modern global navigation satellite system (GNSS) data, and we reduced our surveys to this datum. We surveyed two types of RSL coral proxies (Figs. 1 and 2): The first proxy was fossil microatolls (Porites, Goniastrea, and Heliopora spp.) with annular shapes and flat upper surfaces. Similar to previous studies, we treated the microatoll data as RSL index points that provide upper and lower boundaries (Woodroffe et al., 2012; Kench et al., 2020). The second proxy, fossil in situ corals without microatoll morphology, was interpreted as marine-limiting data indicating that RSL was higher during their growth (Yokoyama and Esat, 2015). The indicative range (IR) and reference water level (RWL) of each index point were determined based on geomorphic observations regarding the possibility of microatoll survival above open-water low-tide levels (i.e., ponding) and the vertical range of living specimens, largely around low tide levels (−0.74 to –0.21 m MSL, n = 43).

Elevations of RSL proxies (10 RSL index points, four marine-limiting data points) and the reef flat geomorphology were measured with differential GNSS. We applied an IR of 0.53 m and a RWL of −0.48 m for nonponded microatolls. For specimens where ponding during their life time could not be excluded, we applied a conservative IR of 1.65 m and RWL of −0.29 m, based on the vertical difference between modeled lowest astronomical tide and mean high water neaps (Egbert and Erofeeva, 2002). Three nonponded fossil microatolls with mostly pristine surface morphologies were slabbed using hand saws and X-rayed, and their internal structures were investigated as RSL proxy indicators following Meltzner and Woodroffe (2015).

Accelerator mass spectrometry radiocarbon ages were determined from aragonitic (X-ray diffraction >93% aragonite) samples. Conventional radiocarbon ages were calibrated to thousands of years B.P. (ka) with OxCal (Ramsey, 2008) using the marine calibration data set Marine20 (Heaton et al., 2020) and a Delta-R value of −64 ± 70 ¹⁴C yr (Southon et al., 2002). Weighted averages of ¹⁴C dates from subsamples of microatoll slab growth bands were determined using the D(·) function in OxCal. We employed a hierarchical statistical model to the RSL proxy indicator data wherein the timing was analyzed through floating chronologies (Ashe et al., 2019).

The glacial rebound solutions used here were developed at the Australian National University (ANU; Nakada and Lambeck, 1988; Lambeck, 1993; Lambeck et al., 2014, 2017). The applied isostatic model focuses on regional solutions for both the rheology and load parameters following an iterative process where resulting viscosity parameters become effective parameters that reflect the response of different parts of the mantle and that are not directly comparable. For this study, we used a three-layer far-field continental margin Earth model based on inversions of a global suite of far-field RSL data reaching back to marine isotope stage 5 (Lambeck et al., 2014). Resulting viscosity parameters suggested an effective elastic lithospheric thickness that ranged between 45 and 80 km, an upper-mantle viscosity that is well constrained to (1.5–2.5) × 10²⁰ Pa·s, and a lower-mantle viscosity for which the upper limit is rather loosely constrained to (1.7 to >20) × 10²² Pa·s. Full details are available in the Supplemental Material.

Upper surfaces of fossil microatolls occupy an elevation range between 0.46 and 1.49 m MSL, whereas fossil corals (marine limiting) have lower elevations from −0.08 to 0.59 m MSL (Table 1). Dated samples provided calibrated age ranges between 6.9 and 6.3 ka. The two oldest, potentially ponded microatolls (PPJ_FMA7 and PPJ_FMA8) are located in a channel behind the rubble rampart and indicate that RSL first reached its present position before ca. 7 ka. Synthesized with recently standardized RSL data from nearby Teluk Awur (Azmy et al., 2010; Mann et al., 2019), the complete record suggests that RSL was ~1.3 ± 0.7 m above present between 6.9 and 5.3 ka (Fig. 3).

This is consistent with predicted RSL because RSL would have risen rapidly during the late deglacial period as a result of the increase in ocean volume. A subsequent change in the rate of ice-volume equivalent sea-level (ESL) rise, coinciding with a period of slowly falling predicted RSL, was caused by the assumed termination of LIS deglaciation at 6.8 ka by the model and reflects the relaxation of the mantle below the ocean crust and the elastic lithosphere response. A period of ~2000 yr followed, during which the ESL contribution from Antarctica of ~2.5 m and loading of the ocean basin adjacent to the study site were nearly balanced.

Significantly, the temporally higher-resolution RSL proxy indicators from coral slabs revealed a RSL trend that would not have been captured if only the less precise fossil microatoll surfaces had been analyzed. Statistically modeled RSL based on age-elevation information from the youngest and oldest growth bands showed a significant and continuous rise by ~1.3 m between 6.7 and 6.4 ka, with maximum significant rates of 6.2 ± 5.0 mm/yr (2σ) between 6.54 and 6.53 ka (Figs. 2 and 3; Table S2).

The occurrence of microatoll RSL index points above coral marine-limiting data suggests a consistent reconstruction presented in this work (Fig. 3). Agreement between the RSL index points and the ANU glacial rebound solution has important implications regarding the potential of predictive isostatic models as tools with which to constrain future RSL variability. First, it is critical to acknowledge the site-specific geological environment that controls the regional isostatic response. Contrary to oceanic islands, for example, at far-field continental margins such as the present study area, the water load is not distributed symmetrically around the site. Accordingly, the isostatic response is a mixture of the mantle response to the distant ice unloading and the ocean response to the water load (Lambeck et al., 2014). Furthermore, the good fit between geological proxies and predictions shows that relatively simple, three-layer mantle solutions as solid-earth model components provide a predictive capability that is consistent with the accuracies of most RSL observations in southeast Asia (Mann et al., 2019).

Proxy indicator data provide a Holocene RSL record on time scales <200 yr that is instructive for resolving ice dynamics. Hierarchical, statistical modeling of the data suggests a sustained RSL rise by ~1.3 m between 6.7 and 6.4 ka, which falls within the 2σ RSL envelope obtained from the individual microatolls that provided time-averaged information. The coherence between the two independent reconstructions and the nearly identical rates of modeled and predicted late-deglacial RSL rise indicate that the most likely explanation for the rapid rise is the final meltwater input from the LIS. Of note, the modeled RSL rise here culminated at 6.4 ± 0.1 ka and thus postdates assumed LIS deglaciation in the applied ice model by 400 yr. While a calibration of this ice model is straightforward due to underlying temporal resolution (500 yr bins), the reconstructed timing of LIS disappearance is significantly younger than that in the ice models proposed by Peltier (2004) and Peltier et al. (2015).

Currently, the most precise chronostratigraphic constraints place ultimate LIS melting between 7.1 and 6.3 ka (Carlson et al., 2008; Ullman et al., 2016). Our interpretation is consistent with these constraints, yet with improved temporal precision, and we suggest that the difference compared to Peltier’s deglaciation chronology is a result of the underlying minimum-limiting radiocarbon data in their ice models.

Interestingly, the peak RSL highstand at 6.4 ± 0.1 ka is ~0.5 m above the predicted median highstand at the end of significant Northern Hemisphere deglaciation (Fig. 3). However, as the reconstructed highstand lies within the confidence limits of the predictions, the difference may be resolved through another combination of Earth-model viscosity parameters. Alternatively, the difference could indicate a larger Antarctic contribution to ESL than is assumed in the ice model. Additional, accurate RSL data across the basin margins, from offshore islands and deeply indented bays and estuaries, are important for analyzing spatial differences between sites in Indonesia to further separate lower- and upper-mantle responses and provide corrective terms to the global ice-volume function (Nakada and Lambeck, 1988).

Financial support for H. Westphal came from the VW foundation through the funding line “Schlüsselthemen,” and for T. Mann from the German Research Foundation (MA 6967/2-1). We thank the Zentrum für moderne Diagnostik (ZEMODI, Bremen) for the X-ray images, Bayu Triyogo Widyantoro and Badan Informasi Geospasial (Indonesia) for providing tide gauge data from Jepara, and RISTEK for providing a research permit (no. 126/SIP/E5/Dit.KI/V/2016). Comments by three reviewers improved the paper.