# High interstadial sea levels over the past 420ka from the Huon Peninsula, Papua New Guinea

Gino de Gelder, Laurent Husson, Anne-Morwenn Pastier, David

Fernández-Blanco, Tamara Pico, Denovan Chauveau, Christine Authemayou,

Kevin Pedoja

► To cite this version:

Gino de Gelder, Laurent Husson, Anne-Morwenn Pastier, David Fernández-Blanco, Tamara Pico, et al.. High interstadial sea levels over the past 420ka from the Huon Peninsula, Papua New Guinea. Communications Earth & Environment, 2022, 3, 10.1038/s43247-022-00583-7. hal-04731748

HAL Id: hal-04731748

https://hal.science/hal-04731748v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# High interstadial sea levels over the past 420ka from the Huon Peninsula, Papua New Guinea

Gino de Gelder ID <fnref n="1" />,²✉, Laurent Husson ID <fnref n="1" />, Anne-Morwenn Pastier³, David Fernández-Blanco ID ⁴, Tamara Pico ID ⁵, Denovan Chauveau ID <fnref n="6" />,⁷, Christine Authemayou<fnref n="6" /> & Kevin Pedoja ID ⁸

The history of sea level across the Quaternary is essential for assessing past and future climate. Global sea-level reconstructions are typically derived from oxygen isotope curves, but require calibration with geological constraints that are scarce prior to the last glacial cycle (>130 thousand years ago). Here we show that the coral reef terrace sequence at the Huon Peninsula (Papua New Guinea) provides such constraints up to ~420 thousand years ago, through a geometric analysis of high-resolution topographic data. We derive a northward tectonic tilt as regional deformation pattern, and estimate relative sea level for 31 Quaternary periods, including several periods for which no relative sea level data exists elsewhere. Supported by numerical reef models, these estimates suggest that oxygen isotope-based global mean sea-level curves systematically underestimate interstadial sea-level elevations, by up to ~20 m. Compared to those curves, our results imply a stronger degree of non-linearity between ice-sheet volumes and global temperatures within Quaternary glacial cycles.

<footnote n="1">¹ISTerre, CNRS, IRD, Université Grenoble-Alpes, St. Martin d'Hères, France. ²Res. Group of Paleoclimate & Paleoenvironment, Res. Centr. for Climate and Atmosphere, Res. Org. of Earth Sciences and Maritime, National Research and Innovation Agency, Bandung, Republic of Indonesia. ³GeoForschungsZentrum, Potsdam, Germany. ⁴Barcelona Center for Subsurface Imaging, Marine Sciences Institute (ICM-CSIC), Barcelona, Spain. ⁵UC Santa Cruz, Santa Cruz, USA.</footnote>

<footnote n="6">⁶UMR 6538 Geo-Ocean, Univ. Brest, CNRS, Ifremer, F-29280 Plouzané, France. ⁷Department of Environmental Sciences, Informatics and Statistics, Ca' Foscari University of Venice, Venice, Italy. ⁸Université de Caen, Caen, France. ✉email: gino.de-gelder@univ-grenoble-alpes.fr</footnote>
Linking past sea-level changes to climatic variations and ice-sheet configurations is fundamental to predicting future climate change scenarios¹ and understanding large-scale geodynamics². The Quaternary represents a key geological archive, as it is the most recent period recording the full spectrum of glacial-interglacial sea-level changes³. Continuous records of past global mean sea level (GMSL) are usually derived from oxygen isotope ratios (δ¹⁸O) in marine sediment cores, as the seawater ratio between the lighter (¹⁶O) and heavier (¹⁸O) isotopes relies on the amount of ice stored within continental ice sheets. Since temperatures and local hydrological/chemical conditions also affect δ¹⁸O ratios, many calibration techniques have been proposed for the derivation of GMSL curves from δ¹⁸O-records, which has led to a broad range of GMSL estimates⁴,⁵. Commonly, to verify or adjust δ¹⁸O-derived sea-level curves, they are compared with reliable relative sea-level (RSL) geological indicators, corrected for Glacial Isostatic Adjustment (GIA). However, these data are sparse⁶,⁷, rely on uncertain corrections for GIA or other processes of vertical displacement, and are mostly limited to interglacial periods⁸ with relatively high sea-level stands - like Marine Isotope Stages (MIS) 5e (~125 ka), 7e (~240 ka), 9e (~320 ka) and 11c (~410 ka)—and interstadial secondary sea-level peaks of the last glacial cycle, like MIS 3 (~45 ka), 5a (~80 ka) and 5c (~100 ka)⁹. Geological RSL indicators of lowstands—the troughs between peaks—and interstadials older than 130 ka are largely unknown or have error bars of 10s of m (Supplementary Fig. 1). This lack of reliable RSL data, combined with uncertainties in δ¹⁸O-based GMSL curves, hinder sea-level, ice-sheet and climatic reconstructions over large portions of the Quaternary.

Here we provide the first geological RSL estimates for several of such periods by re-visiting the canonical coral reef terrace sequence of the Huon Peninsula (Fig. 1; Papua New Guinea). At this location, high uplift rates varying along the coast from ~0.5 mm yr⁻¹ in the NW to 3.5 mm yr⁻¹ in the SE, have led to the preservation of an exceptionally complete sea-level record¹⁰. Whereas previous studies focused on the last glacial cycle terraces, which are easier to access and better preserved, we used high-resolution digital topography to reassess the whole terrace sequence, ~40 km wide and ~1000 m high, with a geometry-based approach. By re-evaluating the tectonic deformation we provide the most complete geological RSL record for the ~420–125 ka period (Supplementary Fig. 1), while improving elevation and age estimates of previously documented sea-level highstands at Huon (~125–0 ka¹¹).

Our data provides a foundation for improved constraints on Quaternary sea-level oscillations and ice-sheets, from which we conclude that δ¹⁸O-derived GMSL-curves systematically underestimate GMSL. Based on this finding, we discuss broader implications for GMSL reconstructions and paleoclimate.

## Results

**New interpretation of Huon coral reef terrace deformation.**
Located on the northern margin of the Huon-Finisterre Ranges on the South Bismarck Plate (Fig. 1), the Huon Peninsula coral reef terraces have provided the basis for fundamental studies on Quaternary sea-level oscillations¹⁰,¹²⁻¹⁶. To interpret RSL in tectonically active regions like the Huon Peninsula, it is crucial to constrain uplift rates and their spatiotemporal variability. To estimate uplift rates, the terrace formed during the MIS 5e peak (~125 ka) is generally used as a benchmark, since the global mean sea-level elevation of MIS 5e is relatively well constrained⁶,¹⁷,¹⁸. At Huon, pioneering studies combined field measurements and observations of the MIS 5e and other terraces to construct detailed uplift rate contours (blue dashed contours in Fig. 1¹⁰),

whereas later studies simply used an uplift rate calculated from the nearest MIS 5e terrace¹⁵,¹⁶,¹⁹. These long-standing interpretations roughly translate as a NW-directed tilt, but this assumption has never been thoroughly verified with a large-scale geometrical analysis. A different tilt direction would imply a different uplift rate correction, and would thus affect sea-level estimates by several meters.

Quantifying km-scale deformation based on field research is complicated, given the lack of large-scale overview at any given location. We better approximate the deformation pattern, and thus sea-level estimations, using a 2 m-resolution Digital Surface Model (DSM) derived from Pleiades satellite images (see Methods) and a 12-m resolution TanDEM-X DSM, which we used to assess the large-scale terrace geometry. Comparison with previously measured profiles²⁰,²¹ shows that our DSM gives similar elevations within a few m (Supplementary Fig. 2). Although field-based estimates of coral reef terraces would typically give smaller error margins for individual transects, the main advantage of using a DSM is that it allows for a continuous evaluation with practically unlimited topographic profiles. This method averages out local peculiarities and avoids the corresponding bias that discrete profiles would give, providing a more objective assessment of terrace elevations. As a qualitative assessment, we used hundreds of parallel topographic swaths of narrow width stacked together perpendicular to their strike²² (hereafter termed stacked swaths²³). Distortion of the terraces would depend on the viewing angle with respect to the dip of the terraces; minimizing the distortion of the terrace sequence indicates a N-ward directed tilt (Supplementary Fig. 3). As a quantitative assessment of this tilt, we systematically measured the elevation of the two most continuous terraces in the sequence (Supplementary Figs. 4 and 5), previously assigned to the MIS 5a and 5e highstands¹⁹. We specifically tracked the shoreline angle elevation, at the intersection between a terrace flat and its corresponding paleo-seacliff landwards (Fig. 2a), as a morphological approximation of paleo-RSL²⁴. Using least-square functions for linear surface fitting, the ~350 measured shoreline angles (Supplementary Data 1) confirm that both terraces are almost exactly N-ward tilted (dip directions N002E for MIS 5a and N359E for MIS 5e; Fig. 2; Supplementary Fig. 5). This tilt direction, ~45 degrees clockwise from the previously assumed NW-directed tilt, better fits the regional tectonic context, given (i) the N-S compression indicated by focal mechanisms of the region's largest earthquakes, (ii) the E-W orientation of the Wonga Thrust immediately S of the Huon Peninsula, and (iii) overall NNE-SSW convergence from GPS constraints²⁵ (Fig. 1). Theoretically, microplate rotation of ~8° per Ma²⁶ could have a minor effect on the deformation pattern through time, but as the deformation pattern does not visibly differ comparing the lower terraces from the upper terraces (Supplementary Fig. 3), we assume that this effect is negligible in the following analysis.

**Inferred relative sea levels.** This improved tectonic interpretation implies different uplift-rate corrections (examples in Supplementary Fig. 2) and, more importantly, different RSL estimates with respect to previous work¹¹. Specifically, compared to a NW-tilt a N-tilt should increase the elevations of RSL estimates younger than MIS 5e, and decrease the elevations older than MIS 5e. We calculate RSL directly from the terrace geometry, taking advantage of the linearity of the shoreline angles. If the deformation pattern and uplift rates have remained constant in time, the slopes of the terraces in across-strike W-ward view (as in Fig. 3a, b) should be proportional to the terrace ages, and extrapolation to a 0 mm·yr⁻¹ uplift rate should give the RSL elevation during terrace formation. We performed this calculation
<figure><img src="image_3.png" /></figure>

**Fig. 1** Uplift rates, dated samples and tectonic context of the Huon coral reef terrace sequence. **a** Slope map of coral reef terrace sequence overlain with colored elevations, showing a comparison of the uplift rates previously proposed¹⁰ and this study. We also show the location of dated samples, and shoreline angles for the upper MIS 5e (~127 ka) terrace. Names in italic refer to the locations of commonly studied sections, of which 3 are modeled in Fig. 4. This map, as well as the ones in **b**, **c**, was made using MAPublisher version 9.8 (https://www.avenza.com/help/mapublisher/9.8/). Basemaps were created with Global Mapper version 15 (https://global-mapper.software.informer.com/15.1/). **b** Regional tectonic setting, with earthquake focal mechanisms of the 1967 and 1992 seismic sequences⁸⁷ ranging from $m_b = 5.5$ (smallest) to 6.2 (largest). GPS arrows²⁵ (white) are drawn with respect to the New Guinea Highlands, as modeled from 66 GPS velocities and 333 earthquake slip vectors. **c** Detailed DSM slope map of the Sialum area, with mapped paleoshorelines for the coral reef terraces. **d** Ages of dated samples within a recent coral sample compilation⁷. All samples in blue, and in orange the samples after removing duplicates, samples out of the study area and U/Th open system samples.

(see Methods) for 31 terraces (T1-T31; Fig. 3b), which we traced using the Pleiades DSM, and constitutes a major expansion from the ~20 previously recognized terraces¹⁰,¹¹. Instead of using a nomenclature based on reef structure¹⁰, on which we have no direct observations, we opted for a simpler designation of terrace numbers (T1-T31; Fig. 3; Supplementary Table 1). To calculate the regional uplift rate pattern we used an age of 127 ± 2 ka, as previously used for the uppermost MIS 5e Huon terrace¹¹ and RSL elevation of 7 ± 5 m. Within the considerable debate on the elevation and timing of MIS 5e sea-level at different locations¹⁷,²⁷⁻³⁰ we find the range of 2 m to 12 m a reasonable compromise between the lower end of early MIS 5e GMSL estimates²⁹,³¹, if GIA effects are negligible at Huon, and the upper end of early MIS 5e GMSL estimates¹⁷,³², if GIA corrections at Huon are on the order of +4 m²⁷. Some RSL lowstand elevations were previously calculated from the Tewai section (Fig. 1¹⁴,³³), which we updated using our new uplift rate estimates.

New RSL estimates for the last glacial cycle, shown in Fig. 3c, are generally higher in elevation than previous Huon estimates¹¹. Our estimates are closer to GIA-corrected GMSL estimates for MIS 3, 5a and 5c interstadials, thus largely accounting for

previous discrepancies between Huon and other sites, which show higher interstadial sea level³⁴⁻³⁶. In terms of timing, our calculated RSL ages are similar to previous studies for Huon's well-dated MIS 3 and 5e highstands (Fig. 3c), but slightly younger for MIS 5a and 5c highstands, which have not been well-constrained by U-Th dating (Fig. 1d). Instead of ~84 ± 2 ka (MIS 5a) and ~107 ± 3 ka (MIS 5c) peak ages¹¹, we find ages of ~80 ± 1 ka and ~101 ± 2 ka, more in line with the peak ages of the GMSL-curve shown in Fig. 3c⁴. Also our re-calculated 5b and 5d lowstands are younger and at higher elevation than previously proposed for Huon¹¹,¹⁴, more in line with the displayed GMSL curve⁴ (Fig. 3c). Generally, our proposed terrace ages are in good agreement with dated Huon samples (Supplementary Table 1). We give updated uplift rate estimates for the locations of the dated samples in Supplementary Data 2.

Comparison of our new RSL estimates for the past ~420 ka (Supplementary Table 1) and the GMSL curve of Spratt and Lisiecki⁴, as displayed in Fig. 3d, shows that highstand ages correspond well, and interglacial highstands are consistently at higher elevations than interstadial highstands. Absolute elevations for the interglacial highstands (5e, 7e, 9e, 11c) are similar for both
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Best-fitting tilt direction of Huon coral reef terrace sequence. **a** Field-observations re-drawn from a previous study¹⁰, with a typical coral reef terraces at the Huon Peninsula, illustrating the relationship between shoreline angle and sea-level at the coastline. **b** Example of shoreline angle determination, at the intersection between terrace and paleo-cliff, using a free-cliff analysis⁷². Black line marks the maximum elevation of the swath profile in Fig. 1c. The orange and blue lines represent a minimum and maximum approximation, respectively, to derive the shoreline angle elevation. **c** Histogram of residuals of surface fitting including standard deviation (STD), with distribution not significantly different from Gaussian at 95% confidence as suggested by a X² test. **d** Linear surface fitting of MIS 5e shoreline angle elevations (blue dots), using the average values of the minima and maxima derived as in **b**. Results indicate a simple planar tilt approximately N-ward (N359E).</figcaption></figure>

GMSL and Huon RSL, and we highlight that the similarity of the Huon MIS 11c RSL highstand to independent GMSL estimates³⁷⁻³⁹ lends support to our assumption that the uplift rate has not changed much over time. More specifically, the difference of ~0–15 ka between the ~400 and 415 ka GMSL highstand age estimate of Dutton et al.³², and our ~415 ka RSL estimate (Fig. 3d) would imply a maximum uplift rate change of 5% between 420 and 125 ka compared to 125-0 ka. In contrast to interglacial sea-levels, interstadial sea-level highstands (MIS 3, 5a, 5c, 6, 9a, 11a) for Huon seem to be of consistently higher elevation than δ¹⁸O-derived GMSL estimates, up to ~20 m for the lower interstadial peaks (Fig. 3d). The offset appears cyclic and scales with absolute sea-level: larger during lower sea levels and

smaller during higher sea levels. We compared our results to an additional 10 δ¹⁸O-derived sea-level curves (Supplementary Fig. 1), obtained with several different methodologies relying on different assumptions, and even though there are exceptions for some interstadials in some curves, the overall trend of interstadial Huon RSL at higher elevations remains.

To compare RSL directly with GMSL, it must be corrected for the local effect of glacial isostatic adjustment (GIA). We used GIA modeling (see Methods) with a global ice history characterized by a GMSL history based on the Spratt and Lisiecki GMSL curve⁴. We found that the maximum GIA correction at this far-field location is less than 7 m, within the uncertainty for most of our RSL estimates (Supplementary Fig. 6). Consequently, our RSL
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Huon relative sea-level estimates. a 500 Stacked swath profiles of TanDEM-X Digital Surface Model looking towards the W with **b** interpretation of coral reef terrace levels. c Zoom-in to last glacial-interglacial cycle, with Huon relative sea-level estimates of this study obtained from coral reef terrace geometry in **b** (Huon RSL - CRT; halo colors as in **b** and perpendicular stacked swath profiles (Supplementary Data 3 and 4; see Methods). We compare these to previous Huon RSL estimates¹¹, global mean sea-level estimates of MIS 3³⁵, MIS 5a³⁶, and MIS 5c³⁶, and a global mean sea-level curve obtained with a Principal Component Analysis⁴ with 97.5% confidence interval (grey band). Lowstand estimates (Huon RSL - Tewai) were re-calculated from other studies¹¹,³³ (see Methods and Supplementary Table 3). Error bars for data points indicate Standard Errors (see Methods). d Same as c but for ~425 ka (halo colors as in **b**, and compared to independently constrained global mean sea-level estimates³² (MIS 5e, 11c). Marine Isotope Stages (MIS; bold) follow a standardized definition⁸⁸.</figcaption></figure>

estimates can be considered to approximate GMSL, and this RSL history constitutes a key dataset, especially for (>130 ka) interstadials, for which we present some of the first sea-level constraints based on geologic indicators.

**Terrace sequence modeling.** Different sea-level histories produce distinct coral reef terrace morphologies, as illustrated by numerical models of coral reef terrace sequences⁴⁰. Models provide a way to precisely incorporate the geomorphic responses of coral reef sequences to sea-level changes, thus allowing for the possibility that rapid sea-level changes and short term fluctuations may not be recorded in fossil reefs⁴¹,⁴². We modelled the Huon coral reef terrace sequence to test (1) whether a full sea-level history with updated uplift rates and highstands adjusted to our RSL data would reproduce the main features of the observed

sequence in a direct model and (2) how much the revised RSL estimates would impact the geomorphic architecture of the sequence. We chose three representative cross-sections (Fig. 4) and used the REEF code⁴²,⁴³. We test the impact of our RSL estimates by considering 3 different sea-level curves (Fig. 4a; Supplementary Fig. 7): (1) a representative GMSL-curve⁴, (2) a modified version of this curve with the whole curve adjusted upwards to match our RSL highstand estimates (3) a similarly modified curve but with the lowstands fixed. To a first order, the model outputs show that our proposed chronology and uplift rates produce reasonable terrace shapes in comparison to the observed sections (Fig. 4b–d). In the model outputs, major highstands (MIS 5e, 7e, 9e, 11c) produce wider terraces than interstadial terraces, a pattern similar to the observed morphology. All dated samples at the Sialum and Kanzarua sections fit with the modeled reef ages. Checking in more detail, two relevant
<figure><img src="image_6.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_9.png" /></figure>

model results stand out: (1) the Sialum model shows how the 119.3 ka sample just below T14 can have an age corresponding to MIS 5e (Fig. 4c), even if the geomorphic structure itself is younger (MIS 5c), and (2) The models also show how multiple terraces can result from a single MIS 5e RSL peak (as in Pastier et al.⁴⁴), and how the same terrace could have been constructed during

multiple times, as observed in the data with multiple clusters of ages between ~90 ka and ~140 ka for the same terrace (Supplementary Table 1;¹⁹,²⁰,⁴⁵).

Considering the misfits between modeled and observed terrace elevations, curves with adjusted highstands better fit the observed terrace morphology, particularly for interstadial terraces T11,
**Fig. 4** Huon terrace sequence modeling. a Sea-level curves used for terrace modeling, including a GMSL curve based on Principal Component Analysis⁴ (black), a modified version of that curve with only the highstands (HS) adjusted (red) and a modified version with both high- and lowstands (HS/LS) adjusted (blue). Corrections were only applied for the shown RSL estimates: given the limited number of peaks in the GMSL-curve we only picked the highest Huon RSL estimate per MIS to calibrate. Error bars for Huon RSL data points indicate Standard Errors (see Methods). b Modeling results for the Gagar Anununai, c Sialum and d Kanzarua sections, using a potential reef growth rate of 10 mm yr⁻¹, erosional potential of 30 mm yr⁻¹ and basal slopes derived from average topography. Misfits are calculated by subtracting the modeled from the observed terrace elevation for the terraces labeled in a. Location of stacked swath profiles, derived from the Pleiades DSM, are given in Fig. 1. The dated samples for the Sialum section were taken within the area of the stacked swath profile⁴⁵, whereas the dated samples for the Kanzarua section¹⁵,⁴⁵,⁵⁸ were taken from a location -2 km S of the stacked swath profile.

T14, T20, T25, T27, T29, and T30 in Fig. 4b, T7, T14, T20, T25, T27, and T29 in Fig. 4c, as well as T7, T11, T14, and T20 in Fig. 4d. Differences between the two highstand-adjusted curves are subtle, and preclude a clear preference. Future work on finer scale may integrate the morphology on m-scale and ages on ka-scale, and incorporate detailed morpho-stratigraphy¹⁶ in order to constrain Huon lowstands using reef modeling.

We also used reef modeling to test two possible limitations that may have affected our RSL estimates. The first possible limitation is that the actual age of construction of a terrace surface may slightly precede sea-level highstands, and change with uplift rate, as was evidence for the Holocene Huon terrace⁴⁶ (Supplementary Fig. 8). The second potential limitation is that a terrace surface might form several meters below the peak in sea-level (SL-Terrace<sub>DIF</sub>), depending on the rate of sea-level changes and on the morphology formed during preceding sea-level oscillations or erosion of a barrier reef¹². We quantified these effects (Supplementary Fig. 8) using the range of uplift rates determined for Huon (0.5–3.5 mm yr⁻¹) in a series of models using the same parameters as for the models in Fig. 4 (see Methods). These reef simulations suggest that peak-delays and non-negligible values of SL-Terrace<sub>DIF</sub> could have lead to misinterpretations in Huon RSL highstand elevation, up to ~10–15 m, and age, up to ~10 ka, for the T31 (MIS 11c) terrace, and that values gradually increase with age. However, as offsets between Huon RSL highstands and δ¹⁸O-derived GMSL estimates appear to be cyclic, rather than increasing gradually (Fig. 3d), neither peak-delays nor important values of SL-Terrace<sub>DIF</sub> can easily explain the highstand discrepancy for interstadial sea-levels.

## Discussion

Collectively, our geometric analysis, and reef modeling results suggest that δ¹⁸O-derived curves systematically underestimate sea level during interstadial periods (Fig. 3; Supplementary Fig. 1). A common argument against using RSL estimates like those in Huon is the possibility of non-constant uplift rates⁴⁷,⁴⁸. While we cannot exclude this possibility, the cyclic pattern of large discrepancies with interstadial highstands and low discrepancies with interglacial highstands (Fig. 3d; Supplementary Fig. 7) is not easily reconciled with a tectonic periodicity. Even though ~0–3 m episodic uplift events at timesteps of 200–1900 years have probably occurred throughout the Late Pleistocene⁴⁹, it seems unlikely that clusters of tectonic uplift events would systematically occur much more frequent during interstadial periods. As such, the simplest explanation for the consistently higher Huon RSL estimates in comparison to δ¹⁸O-derived GMSL estimates is that the latter does not appropriately account for reduced ice-sheet extent of interstadial highstands. Combined geologic/paleogeographic observations and GIA analysis suggest that indeed global ice volumes may be smaller than inferred from δ¹⁸O records during MIS 3³⁵,⁵⁰–⁵². Similar trends may hold for the older interstadials peaks within MIS 6, 9a, 11a. Previous studies highlighted several problems in the conversion from δ¹⁸O to sea-level⁴,⁵³–⁵⁵. Along the lines of those studies, we envision two possibilities, either (1) δ¹⁸O has been incorrectly calibrated to sea-level in previous work,

or (2) δ¹⁸O does not fully grasp short-lived changes in the growth and decay of ice sheets to begin with.

Considering option (1) above, Waelbroeck et al.⁵⁵ proposed that δ¹⁸O should be calibrated differently for deglaciations than for glaciations. They performed a separate calibration with RSL indicators, and this is probably why their interstadial highstand elevations are more similar to the Huon data than most other δ¹⁸O-derived GMSL curves (Supplementary Fig. 1), even though they relied on the previous, lower Huon RSL estimates for the interstadial highstand periods. An updated calibration with a similar technique could reconcile δ¹⁸O with geologic indicators. An alternative explanation for relatively low sea-level estimates in δ¹⁸O-derived curves is the possibility that sea-level during the last glacial maximum (LGM), to which most δ¹⁸O-derived curves are calibrated, may be generally estimated too low⁵⁶.

Considering option (2) above, it is possible that rapid and high-amplitude sea-level oscillations are long enough to imprint uplifting landscapes, but too short to be well recorded in δ¹⁸O records. Such records are often stacked or averaged, ignoring potentially important differences in timing between oceanic basins⁵⁵,⁵⁷, and smoothing out rapid sea-level changes⁴. Rapid (~7 ky periodicity), high-amplitude (10–15 m) oscillations have been proposed for MIS 3 in Huon¹⁶,⁵⁸ in relation to Heinrich events⁵⁹, supported by the large number of Huon MIS 3 terraces and derived RSL estimates (Fig. 3). Heinrich events may also exist for older interstadial periods than MIS 3⁶⁰, although this has never been directly evidenced by geologic data. We note in our analysis that, similar to T2–T8 for MIS 3, T17–19 seem to represent a ~7 ky periodicity for the MIS 6 interval (Fig. 3d). A detailed morpho-stratigraphic analysis of these possible MIS 6 terraces, in combination with dating and modeling¹⁶,⁵⁸ may shed more light on the possibility of older Heinrich events.

To place our results in a broader climatic perspective, we compare Quaternary GMSL and Huon RSL with paleotemperature, CO₂, and northern hemisphere summer insolation estimates (Fig. 5). Although the continental paleotemperature record (Fig. 5a) is a relative record for Antarctica, it should be representative of global mean variations, albeit more sensitive with a factor 1.4⁶¹. After each interglacial peak, temperature and, to a lesser extent, CO₂, tend to decrease rapidly early and slowly later in the cycle, whereas for GMSL, this is opposite (transparent arrows in Fig. 5a). This non-linear relationship between GMSL and paleotemperature has been described before, and can be attributed to (a combination of) delayed glacial isostatic rebound⁶²–⁶⁴, differences in glacial thresholds of the different ice-sheets⁶⁵,⁶⁶, and cooling of the globe as a precondition for ice sheet growth⁶⁷. Whichever of these mechanisms is dominant to the relationship between GMSL and paleotemperature/CO₂ (solid arrows in Fig. 5a), our Huon RSL estimates suggest that the effect of non-linear processes has been more pronounced than previously imagined. Considering the principal drivers behind GMSL changes, Shakun et al.⁶⁷ showed, using a global stack of paired sea surface temperature-planktonic δ¹⁸O records, that rates of GMSL-changes are largely accounted for by the combination of N-hemisphere summer insolation, sea-surface temperatures (SST;
<figure><img src="image_10.png" /><figcaption>**Fig. 5** Comparison of paleo sea-level and paleo climate data. **a** Comparison of GMSL⁴, Huon RSL, temperature⁸⁹ and CO₂⁹⁰. Temperature anomaly and CO₂ record are both derived from Antarctic EPICA Dome C ice cores. Vertical axes are scaled to the maximum and minimum of each record within the last glacial cycle (~130 ka - present). Transparent arrows illustrate approximate trends in temperature and GMSL per glacial cycle, solid arrows illustrate the stronger non-linearity in GMSL as suggested by Huon RSL constraints. Error bars for Huon RSL data points indicate Standard Errors (see Methods). **b** Insolation variations at 65° N on June 21<sup>91</sup> and sea-surface temperature variations from a global stack⁶⁷. **c** Rate of change of the sea-level curves in **a**, calculated over a running 6-kyr window and scaled to a multivariate regression model (as in previous work⁶⁷) based on the summer insolation and global sea-surface temperatures in **b**, and sea-level curves in **a**.</figcaption></figure>

Fig. 5b) and absolute GMSL. We applied a similar multivariate regression (Fig. 5c; Supplementary Fig. 9) to the 3 sea-level curves of Fig. 4a, corroborating that insolation, SST, and GMSL can explain ~70% of the variability, but note that adjusting the highstands based on Huon RSL only alters the outcomes of the regression to a minor extent. We conclude that our RSL data from the Huon Peninsula provides a basis to accurately constrain Quaternary sea-level and its relation with paleoclimate, ice-sheets, and geodynamics, motivating a re-evaluation of interstadial highstands and relations of sea-level and δ¹⁸O.

**Methods**

**Digital surface model (DSM) calculation.** We used tri-stereo 0.5 m-resolution Pleiades satellite images to calculate the topography of the Huon Peninsula terrace sequence. The open-source software Ames Stereopipeline⁶⁸ allowed us to produce ortho-rectified images, point clouds and a 2 m-resolution DSM, to which we applied a ~58.9 m vertical correction to account for the geoid height. We ground-checked absolute elevations of 10 control points along the coastline, that were all between 0 and 1 m elevation. We used a 30 m-resolution ALOS DSM to fill in the gaps that resulted from cloud coverage within the Pleiades DSM.

**Terrace tilt.** We calculated large-scale parallel swath profiles with Topobox 2.⁶⁹ and stacked them orthogonal to their trend using an in-house script. This wide swath²², later termed stacked swaths²³, is a 2D profile of hundreds of narrow-width topographic swath profiles plot together as thin-air lines. Stacked swaths highlight the

coherence of the topography along the point of view, i.e. in depth, and the consistency in slope and overall morphology of the topography perpendicular to the point of view. We find the view of the topography that is most coherent regarding the tilt of the terrace sequence, i.e. that perpendicular to the tilt (Supplementary Fig. 3), using the average topography of swath profiles on the TanDEM-X DSM and several viewing angles. The width of each swath profile is computed dynamically, as the total width of the data in the viewpoint direction divided by the chosen number of swaths. We secure an average swath profile width of ~50 m, regardless of viewing angle, plotting a variable number (between 400 and 700) of swath profiles.

**Terrace deformation.** As terraces are not always flat (e.g. example in Fig. 2b), either because of the way they formed and/or how they were eroded afterwards, we chose to focus our evaluation of the deformation on shoreline angles (Fig. 2). We followed the detailed mapping of the Sialum section in Stein et al.²⁰ and age interpretations of Esat et al.¹⁹ to locate the continuous MIS 5a and 5e terraces on the DSM. From this location, we continued mapping paleoshorelines of these terrace levels NW and SE ward along the coastline. For this, we used a combination of satellite imagery, slope maps and hillshade images of the DSM. We derived shoreline angles by systematically placing ~350 cliff-perpendicular swath profiles along the inner edges of the MIS 5a and 5e terraces, using TerraceM⁷⁰ (to Supplementary Fig. 4). To reduce the influence of gravitational, karstic and/or fluvial erosion we used the maximum topography of 100 m-wide swath profiles (as in Jara-Muñoz et al.⁷¹). Assuming a broad range of 60–90° as paleocliff slope, we used a free-cliff analysis²² to calculate minimum and maximum shoreline angle elevations by picking a most seaward/landward position of the paleo-seacliff, respectively (Fig. 2a; Supplementary Data 1). The shoreline angles plotted in Figs. 1 and 2 are the average values between minimum and maximum calculated positions. We
used a linear surface fitting function to constrain the deformation pattern, assuming a simple, linear tilt, for both the MIS 5a and 5e shoreline angles (Fig. 2; Supplementary Fig. 5). We performed critical χ2 tests to confirm that the residuals follow a Gaussian distribution and the shoreline angle dataset is well described by a linear surface, which was the case for both the MIS 5a and 5e shoreline angles (Fig. 2; Supplementary Fig. 5).

**RSL highstand calculation.** Coral reef terraces in Huon were previously described with Latin numerals and subdivisions¹⁰. Given that we do not have direct observations of the internal reef structure, we opted for a simpler designation of terrace numbers (T1-T31; Fig. 3; Supplementary Table 1). To calculate uplift rates, *U*, we used the shoreline angle elevations of T16 (MIS 5e), as follows:

$$U = \frac{H_{SA} - H_{RSL} + H_{REF}}{T} \qquad (1)$$

where $H_{SA}$ is the shoreline angle elevation, $H_{RSL}$ the relative sea-level elevation with respect to present-day sea-level, $H_{REF}$ the average reference water level for Huon coral reef terraces and *T* the time of terrace formation. To calculate standard errors we used:

$$\sigma_U^2 = U^2 \left( \left( \frac{\sigma_{HSA}^2 + \sigma_{HRL}^2 + \sigma_{HREF}^2 + \sigma_{HSEIS}^2 + \sigma_{HDSM}^2}{(H_{SA} - H_{RSL} + H_{REF})^2} \right) + \left( \frac{\sigma_T^2}{T^2} \right) \right) \quad (2)$$

where $\sigma_{HSA}$, $\sigma_{HRL}$, $\sigma_{HREF}$, $\sigma_{HSEIS}$, $\sigma_{HDSM}$, and $\sigma_T$ are the uncertainties in shoreline angle elevation, RSL, reference water level, seismic cycle influence, DSM elevation and time of terrace formation, respectively. We calculated $\sigma_{HSA}$ by dividing the difference between the maximum and minimum shoreline angle estimates by two. Typical differences between minimum and maximum shoreline angle elevations are on the order of ~2–10 m, with some extremes of a few 10 s of m for the older coral reef terraces (Supplementary Data 3, 4). For $H_{RSL}$ and $\sigma_{HRL}$ we use $7 \pm 5$ m here, as a compromise between reasonable MIS 5e minimum and maximum estimates of 2 m³ and 12 m²⁷, respectively. For $H_{REF}$ and $\sigma_{HREF}$ we use $0.66 \pm 0.33$ m, as calculated with the indicative meaning calculator⁷³ (Supplementary Table 2). For $\sigma_{HSEIS}$ we used the results of Ota and Chappell⁴⁹ on the Holocene seismic cycles, who found that at any given moment the short-term uplift was within ±1.5 m of the long-term expected uplift at a site with long-term averaged ~3 mm yr⁻¹ uplift rate. Assuming that long-term uplift is a cumulative effect of several seismic cycles⁷⁴ and thus the long-term uplift rate is proportional to the seismic cycle deviation, we use $\sigma_{HSEIS} = U/2$ as a rough estimation. For $\sigma_{HDSM}$ we use a value of 1 m, which is a conservative estimate given that other studies using the same calculation method with Pleiades satellite images found DSM standard errors of 0.3–0.5 m⁷⁵⁻⁷⁷. For *T* and $\sigma_T$ we use $127 \pm 2$ ka for MIS 5e, as was used in Lambeck and Chappell¹¹ following assessment of dated samples from both Huon and other sites (see their paper for details).

We express the uplift rate as a linear function of latitude *l* in meters to propose highstand ages and elevations for the non-dated terraces. Following our deformation assessment, we used a weighted fit of the shoreline-angle derived uplift rate estimates in the form $U = A + Bl$. Following Taylor⁷⁸, we introduced the weights $w = 1/\sigma_U^2$, with best estimates for constants *A* and *B* calculated from:

$$A = \frac{\sum wl^2 \sum wU - \sum wl \sum wU}{\Delta} \qquad (3)$$

$$B = \frac{\sum w \sum wU - \sum wl \sum wU}{\Delta} \qquad (4)$$

with

$$\Delta = \sum w \sum wl^2 - (\sum wl)^2 \qquad (5)$$

And the uncertainties in *A* and *B* as:

$$\sigma_A = \sqrt{\frac{\sum wl^2}{\Delta}} \qquad (6)$$

$$\sigma_B = \sqrt{\frac{\sum w}{\Delta}} \qquad (7)$$

We used the calculated relation between *U* and *l* to estimate the latitude (given a strict northward tilt) for which $U = 0$ mm yr⁻¹, and subsequently to estimate RSL highstand elevations and ages for the other coral reef terraces. We quantified shoreline angle elevations for all terraces from E-W oriented, 2 m resolution stacked swath profiles (Supplementary Data 3), with widths of ~700 m corresponding to uplift rate steps of approximately 0.05 mm yr⁻¹ (Supplementary Fig. 4). We then added $H_{REF}$ of 0.66 m to every shoreline angle estimate to obtain past sea-level elevations $SL_p$, with a standard error calculated from $\sigma_{HSA}$, $\sigma_{HREF}$, $\sigma_{HSEIS}$ and $\sigma_{HDSM}$. Using Formulas 3–7, but replacing *U* with $SL_p$, we obtained a linear weighted fit of the corrected shoreline angles for every terrace (Supplementary Data 4). We estimate the age of each terrace, using pure proportionality:

$$T_{TER} = T_{T16} \cdot \frac{B_{TER}}{B_{T16}} \qquad (8)$$

with standard error

$$\sigma_{T_{TER}} = T_{TER} \cdot \sqrt{\left(\frac{\sigma_{T_{T16}}}{T_{T16}}\right)^2 + \left(\frac{\sigma_{B_{TER}}}{B_{TER}}\right)^2 + \left(\frac{\sigma_{B_{T16}}}{B_{T16}}\right)^2} \qquad (9)$$

Where $T_{T16}$, $B_{T16}$ and its corresponding standard errors are age and slope, respectively, from the T16 (MIS 5e) terrace. To estimate the RSL elevation for each terrace, we used

$$H_{RSL} = A_{TER} + (B_{TER} \cdot l_{u=0}) \qquad (10)$$

with standard error

$$\sigma_{H_{RSL}} = \sqrt{\sigma_{A_{TER}}^2 + \left( B_{TER} \cdot l_{u=0} \cdot \sqrt{\left( \frac{\sigma_{B_{TER}}}{B_{TER}} \right)^2 + \left( \frac{\sigma_{l_{u=0}}}{l_{u=0}} \right)^2} \right)^2} \qquad (11)$$

In which $A_{TER}$ and $B_{TER}$ are constants to fit each terrace. Finally, to re-calculate the lowstand elevations given in previous work¹¹,¹⁴,³³ we considered the elevations of delta deposits in the Tewai section (Fig. 1), but we used our new uplift rate estimates instead of previously used uplift rates of 3.3–3.5 mm yr⁻¹ (Supplementary Table 3). We use the same age and elevation uncertainties as previously given to these lowstands, and as previously done, place our best guess for the age in the middle between the preceding and following highstands.

**GIA calculation.** To calculate relative sea-level change at the Huon Peninsula, we used a gravitationally self-consistent glacial isostatic adjustment model. Our calculations are based on the theory and pseudo-spectral algorithm described by Kendall et al.⁷⁹ with a spherical harmonic truncation at degree and order 512. These calculations solve the full sea-level equation, including the impact of load-induced Earth rotation changes on sea level⁸⁰,⁸¹, evolving shorelines and the migration of grounded portions of marine-based ice sheets⁷⁹,⁸²⁻⁸⁴.

Our numerical predictions required models for Earth’s viscoelastic structure and the history of global ice cover. We adopted the viscosity profile VM²⁸⁵, and tested the sensitivity to this choice by running two simulations using earth models characterised by (1) a lithospheric thickness of 71 km, an upper and lower mantle viscosity of 3 × 10²⁰ Pa s, and 5 × 10²¹ Pa s, respectively, and (2) a lithospheric thickness of 48 km, an upper and lower mantle viscosity of 5 × 10²⁰ Pa s, and 15 × 10²¹ Pa s, respectively. These earth structure parameters are in agreement with the range estimated in Lambeck and Chappell¹¹. We created a global ice sheet history associated with the GMSL history inferred in Spratt and Lisiecki⁴. We modified the GMSL curves to include no excess melting during interglacials. For GMSL values greater than zero, we assigned the present-day ice sheet geometry with a GMSL value of zero. We constructed this ice history by assuming that ice geometry in the pre-LGM period is identical to the time in the post-LGM period with the same GMSL value, using the deglacial ice sheet reconstruction in ICE-5G⁸⁵. Our glacial isostatic adjustment calculations were performed for ice sheet loading changes at increments of 2 ky from 430 ka to 32 ka, increments of 1 ky from 32 ka to 17 ka, and increments of 0.5 ka from 17 ka to 0 ka (similar time steps to ICE-5G).

**Terrace sequence modeling.** We opted for the GMSL curve of Spratt and Lisiecki⁴ without the GIA-correction, as the differences with the corrected curve are very small (<5 m, Supplementary Fig. 6), and this way the effect of highstand-adjustment is slightly more pronounced and easier to observe. We used shape-preserving spline interpolations to calculate the adjusted sea-level curves (Supplementary Fig. 7), modified from the GMSL-curve of Spratt and Lisiecki⁴. For the highstand-adjusted curve we used the difference between our calculated RSL highstands and the GMSL highstands, and values of 0-m for lowstands, whereas for the high- and lowstand adjusted curves we left out the 0-m values for lowstands (Supplementary Fig. 7). For the MIS peaks characterised by multiple RSL estimates, we picked the highest RSL estimate to adjust to the GMSL-curve of Spratt and Lisiecki⁴. For the reef modeling, we used a potential reef growth rate of 10 mm yr⁻¹, as was used by Chappell¹⁶ based on U-series dating of the Late-Glacial reef⁸⁶. We used a basal slope similar to the average slope along the modelled profiles, and erosional potential of 30 mm yr⁻¹ (see full description of the parameterization in Pastier et al.⁴²).

We did a series of modeling tests using the highstand adjusted sea-level curve to check the possible effects of sea-level peak-delays and coral reef terraces forming several meters below sea-level (SL-Terrace<sub>DIF</sub>; Supplementary Fig. 8). As for Fig. 4, we modelled terrace sequences for an uplift rate range of 0.5–3.5 mm yr⁻¹, potential reef growth rate of 10 mm yr⁻¹, erosional potential of 30 mm yr⁻¹ and basal slopes similar to Huon stacked swath profiles (Supplementary Data 3). To constrain the peak-delays, we checked the ages at which the shoreline angles of the different coral reef terraces formed, and the difference between the peak in sea-level of the corresponding highstand (Supplementary Fig. 8a, c). To constrain SL-Terrace<sub>DIF</sub>, we checked the difference between the theoretical elevation of the coral reef terrace and the elevation at which the terrace formed (Supplementary Fig. 8b, d). Following the same approach as in the previous section (RSL highstand calculation; Fig. 2), we then calculated RSL highstands from the morphology of the coral reef terrace sequence, to compare with the RSL curve used as model input. This allowed us to have a rough
estimate of how much peak-delays and SL-Terrace<sub>DIF</sub> could have influenced our Huon RSL estimates (Supplementary Fig. 8e).

## Data availability

The Pleiades satellite imagery was obtained through the DINAMIS program of the Centre National d'Etudes Spatiales (CNES, France) under an academic license and is not for open distribution. We do however share the Digital Surface Model that was derived from the satellite imagery (Supplementary Data 5), in GeoTiff format with 2-m horizontal resolution, which is enough to reproduce the results of this study. We note that the original imagery could be requested through dinamis@cnes.fr (please quote this paper). The 12-m resolution TanDEM-X Digital Surface Model that is used in this study was also obtained under an academic license and not for open distribution. It is not necessary to possess this data in order to reproduce the results of this study: this data was only used for aesthetic purposes in some of the maps, and Supplementary Fig. 3. To freely obtain such data, we note that calls for scientific projects are regularly launched through the TanDEM-X Science website (https://tandemx-science.dlr.de/).

## Code availability

The code for reef modeling is available through https://github.com/Anne-Morwenn/REEF. The TerraceM-2 code for shoreline angle calculation is available through https://doi.org/10.3389/feart.2019.00255.

Received: 28 August 2021; Accepted: 12 October 2022; Published online: 29 October 2022

## References

1. DeConto, R. M. & Pollard, D. Contribution of Antarctica to past and future sea-level rise. *Nature* **531**, 591–597 (2016).

2. Conrad, C. P. & Husson, L. Influence of dynamic topography on sea level and its rate of change. *Lithosphere* (2009).

3. Murray-Wallace, C. V. & Woodroffe, C. D. *Quaternary Sea-Level Changes: A Global Perspective*. (Cambridge University Press, 2014).

4. Spratt, R. M. & Lisiecki, L. E. A Late Pleistocene sea level stack. *Clim. Past* **12**, 1079 (2016).

5. de Gelder, G. et al. How do sea-level curves influence modeled marine terrace sequences? *Quat. Sci. Rev.* **229**, 106132 (2020).

6. Pedoja, K. et al. Coastal staircase sequences reflecting sea-level oscillations and tectonic uplift during the Quaternary and Neogene. *Earth Sci. Rev.* **132**, 13–38 (2014).

7. Hibbert, F. D. et al. Coral indicators of past sea-level change: A global repository of U-series dated benchmarks. *Quat. Sci. Rev.* **145**, 1–56 (2016).

8. Past Interglacials Working Group of PAGES. Interglacials of the last 800,000 years. *Rev. Geophys.* **54**, 162–219 (2016).

9. Thompson, S. B. & Creveling, J. R. A global database of marine isotope substage 5a and 5c marine terraces and paleoshoreline indicators. *Earth Syst. Sci. Data* **13**, 3467–3490 (2021).

10. Chappell, J. Geology of Coral Terraces, Huon Peninsula, New Guinea: A Study of Quaternary Tectonic Movements and Sea-Level. *Changes. GSA Bulletin* **85**, 553–570 (1974).

11. Lambeck, K. & Chappell, J. Sea level change through the last glacial cycle. *Science* **292**, 679–686 (2001).

12. Veeh, H. H. & Chappell, J. Astronomical theory of climatic change: support from new Guinea. *Science* **167**, 862–865 (1970).

13. Bloom, A. L., Broecker, W. S., Chappell, J. M. A., Matthews, R. K. & Mesolella, K. J. Quaternary sea level fluctuations on a tectonic coast: New 230Th234U dates from the Huon Peninsula, New Guinea. *Quat. Res.* **4**, 185–205 (1974).

14. Chappell, J. & Shackleton, N. J. Oxygen isotopes and sea level. *Nature* **324**, 137–140 (1986).

15. Chappell, J. et al. Reconciliation of late Quaternary sea levels derived from coral terraces at Huon Peninsula with deep sea oxygen isotope records. *Earth Planet. Sci. Lett.* **141**, 227–236 (1996).

16. Chappell, J. Sea level changes forced ice breakouts in the Last Glacial cycle: new results from coral terraces. *Quat. Sci. Rev.* **21**, 1229–1240 (2002).

17. Kopp, R. E., Simons, F. J., Mitrovica, J. X., Maloof, A. C. & Oppenheimer, M. Probabilistic assessment of sea level during the last interglacial stage. *Nature* **462**, 863–867 (2009).

18. Rovere, A. et al. The analysis of Last Interglacial (MIS 5e) relative sea-level indicators: Reconstructing sea-level in a warmer world. *Earth-Sci. Rev.* **159**, 404–427 (2016).

19. Esat, T. M., McCulloch, M. T., Chappell, J., Pillans, B. & Omura, A. Rapid fluctuations in sea level recorded at huon peninsula during the penultimate deglaciation. *Science* **283**, 197–201 (1999).

20. Stein, M. et al. TIMS U-series dating and stable isotopes of the last interglacial event in Papua New Guinea. *Geochim. Cosmochim. Acta.* **57**, 2541–2554 (1993).

21. Chappell, J., Ota, Y. & Berryman, K. Late quaternary coseismic uplift history of Huon Peninsula, Papua New Guinea. *Quat. Sci. Rev.* **15**, 7–22 (1996).

22. Armijo, R., Lacassin, R., Coudurier-Curveur, A. & Carrizo, D. Coupled tectonic evolution of Andean orogeny and global climate. *Earth-Sci. Rev.* **143**, 1–35 (2015).

23. Fernández-Blanco, D., Gelder, G., Lacassin, R. & Armijo, R. Geometry of flexural uplift by continental rifting in Corinth, Greece. *Tectonics* https://doi.org/10.1029/2019TC005685 (2019).

24. Lajoie, K. R. Coastal tectonics. *Active tectonics* (1986).

25. Wallace, L. M. et al. GPS and seismological constraints on active tectonics and arc-continent collision in Papua New Guinea: Implications for mechanics of microplate rotations in a plate boundary zone. *J. Geophys. Res.* **109**, (2004).

26. Tregoning, P. et al. Motion of the south Bismarck plate, Papua New Guinea. *Geophys. Res. Lett.* **26**, 3517–3520 (1999).

27. Creveling, J. R., Mitrovica, J. X., Hay, C. C., Austermann, J. & Kopp, R. E. Revisiting tectonic corrections applied to Pleistocene sea-level highstands. *Quat. Sci. Rev.* **111**, 72–80 (2015).

28. Barlow, N. L. M. et al. Lack of evidence for a substantial sea-level fluctuation within the Last Interglacial. *Nat. Geosci.* **11**, 627–634 (2018).

29. Dyer, B. et al. Sea-level trends across The Bahamas constrain peak last interglacial ice melt. *Proc. Natl. Acad. Sci. USA.* **118**, (2021).

30. Rohling, E. J. et al. Asynchronous Antarctic and Greenland ice-volume contributions to the last interglacial sea-level highstand. *Nat. Commun.* **10**, 5040 (2019).

31. Hearty, P. J., Hollin, J. T., Neumann, A. C., O'Leary, M. J. & McCulloch, M. Global sea-level fluctuations during the Last Interglaciation (MIS 5e). *Quat. Sci. Rev.* **26**, 2090–2112 (2007).

32. Dutton, A. et al. Sea-level rise due to polar ice-sheet mass loss during past warm periods. *Science* **349**, aaa4019 (2015).

33. Chappell, J. A revised sea level record for the last 300,000 years from Papua New Guinea. *Search.* **14**, 99–101 (1983).

34. Cabioch, G. & Ayliffe, L. K. Raised Coral Terraces at Malakula, Vanuatu, Southwest Pacific, Indicate High Sea Level During Marine Isotope Stage 3. *Quat. Res.* **56**, 357–365 (2001).

35. Pico, T., Mitrovica, J. X., Ferrier, K. L. & Braun, J. Global ice volume during MIS 3 inferred from a sea-level analysis of sedimentary core records in the Yellow River Delta. *Quat. Sci. Rev.* **152**, 72–79 (2016).

36. Creveling, J. R., Mitrovica, J. X., Clark, P. U., Waelbroeck, C. & Pico, T. Predicted bounds on peak global mean sea level during marine isotope stages 5a and 5c. *Quat. Sci. Rev.* **163**, 193–208 (2017).

37. Raymo, M. E. & Mitrovica, J. X. Collapse of polar ice sheets during the stage 11 interglacial. *Nature* **483**, 453–456 (2012).

38. Roberts, D. L., Karkanas, P., Jacobs, Z., Marean, C. W. & Roberts, R. G. Melting ice sheets 400,000 yr ago raised sea level by 13m: past analogue for future trends. *Earth Planet. Sci. Lett.* **357–358**, 226–237 (2012).

39. Chen, F. et al. Refining estimates of polar ice volumes during the MIS11 interglacial using sea level records from South Africa. *J. Clim.* **27**, 8740–8746 (2014).

40. Toomey, M., Ashton, A. D. & Taylor Perron, J. Profiles of ocean island coral reefs controlled by sea-level history and carbonate accumulation rates. *Geology* **41**, 731–734 (2013).

41. Camoin, G. F. & Webster, J. M. Coral reef response to Quaternary sea-level and environmental changes: state of the science. *Sedimentology* **62**, 401–428 (2015).

42. Pastier, A.-M. et al. Genesis and architecture of sequences of quaternary coral reef terraces: insights from numerical models. *Geochem. Geophys. Geosyst.* **20**, 4248–4272 (2019).

43. Husson, L. et al. Reef carbonate productivity during quaternary sea level oscillations. *Geochem. Geophys. Geosyst.* **19**, 1148–1164 (2018).

44. Pastier, A.-M., Malatesta, L., Huppert, K. & Chauveau, D. Towards a dynamic approach of sequences of coral reef terraces. *10.5194/egusphere-egu21-15976 (EGU General Assembly, 2021).*

45. Cutler, K. B. et al. Rapid sea-level fall and deep-ocean temperature change since the last interglacial period. *Earth Planet. Sci. Lett.* **206**, 253–271 (2003).

46. Ota, Y. & Chappell, J. Holocene sea-level rise and coral reef growth on a tectonically rising coast, Huon Peninsula, Papua New Guinea. *Quat. Int.* **55**, 51–59 (1999).

47. Saillard, M. et al. Andean coastal uplift and active tectonics in southern Peru: 10Be surface exposure dating of differentially uplifted marine terrace sequences (San Juan de Marcona, ~15.4°S). *Geomorphology* **128**, 178–190 (2011).

48. Walker, R. T. et al. Rapid mantle-driven uplift along the Angolan margin in the late Quaternary. *Nat. Geosci.* **9**, 909 (2016).

49. Ota, Y. & Chappell, J. Late Quaternary coseismic uplift events on the Huon Peninsula, Papua New Guinea, deduced from coral terrace data. *J. Geophys. Res.* **101**, 6071–6082 (1996).

50. Batchelor, C. L. et al. The configuration of Northern Hemisphere ice sheets through the Quaternary. *Nat. Commun.* **10**, 3713 (2019).
51. Dalton, A. S. et al. Was the Laurentide Ice Sheet significantly reduced during Marine Isotope Stage 3? *Geology* **47**, 111–114 (2019).

52. Gowan, E. J. et al. A new global ice sheet reconstruction for the past 80 000 years. *Nat. Commun.* **12**, 1199 (2021).

53. Shackleton, N. J. Oxygen isotopes, ice volume and sea level. *Quat. Sci. Rev.* **6**, 183–190 (1987).

54. Adkins, J. F., McIntyre, K. & Schrag, D. P. The Salinity, Temperature, and δ18O of the Glacial Deep Ocean. *Science* **298**, 1769–1773 (2002).

55. Waelbroeck, C. et al. Sea-level and deep water temperature changes derived from benthic foraminifera isotopic records. *Quat. Sci. Rev.* **21**, 295–305 (2002).

56. Simms, A. R., Lisiecki, L., Gebbie, G., Whitehouse, P. L. & Clark, J. F. Balancing the last glacial maximum (LGM) sea-level budget. *Quat. Sci. Rev.* **205**, 143–153 (2019).

57. Skinner, L. C. & Shackleton, N. J. An Atlantic lead over Pacific deep-water change across Termination I: implications for the application of the marine isotope stage stratigraphy. *Quat. Sci. Rev.* **24**, 571–580 (2005).

58. Yokoyama, Y., Esat, T. M. & Lambeck, K. Coupled climate and sea-level changes deduced from Huon Peninsula coral terraces of the last ice age. *Earth Planet. Sci. Lett.* **193**, 579–587 (2001).

59. Heinrich, H. Origin and consequences of cyclic ice rafting in the Northeast Atlantic Ocean during the past 130,000 years. *Quat. Res.* **29**, 142–152 (1988).

60. Hodell, D. A., Channell, J. E. T., Curtis, J. H., Romero, O. E. & Röhl, U. Onset of 'Hudson Strait' Heinrich events in the eastern North Atlantic at the end of the middle Pleistocene transition (~640 ka)? *Paleoceanography* **23**, 4 (2008).

61. Rohling, E. J., Medina-Elizalde, M., Shepherd, J. G., Siddall, M. & Stanford, J. D. Sea surface and high-latitude temperature sensitivity to radiative forcing of climate over several glacial cycles. *J. Clim.* **25**, 1635–1656 (2012).

62. Oerlemans, J. Model experiments on the 100,000-yr glacial cycle. *Nature* **287**, 430–432 (1980).

63. Pollard, D. A simple ice sheet model yields realistic 100 kyr glacial cycles. *Nature* **296**, 334–338 (1982).

64. Abe-Ouchi, A. et al. Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume. *Nature* **500**, 190–193 (2013).

65. de Boer, B., van de Wal, R. S. W., Lourens, L. J. & Bintanja, R. Transient nature of the Earth's climate and the implications for the interpretation of benthic δ18O records. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **335–336**, 4–11 (2012).

66. Gasson, E., Siddall, M. & Lunt, D. J. Exploring uncertainties in the relationship between temperature, ice volume, and sea level over the past 50 million years. *Reviews of Geophysics* **50**, 1 (2012).

67. Shakun, J. D., Lea, D. W., Lisiecki, L. E. & Raymo, M. E. An 800-kyr record of global surface ice δ18O and implications for ice volume-temperature coupling. *Earth Planet. Sci. Lett.* **426**, 58–68 (2015).

68. Beyer, R. A., Alexandrov, O. & McMichael, S. The Ames stereo pipeline: NASA's open source software for deriving and processing terrain data. *Earth Space Sci.* **5**, 537–548 (2018).

69. Schwanghart, W. & Scherler, D. TopoToolbox 2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences. *Earth Surface Dynamics* **2**, 1–7 (2014).

70. Jara-Muñoz, J., Melnick, D., Pedoja, K. & Strecker, M. R. TerraceM-2: A Matlab® interface for mapping and modeling marine and lacustrine terraces. *Front Earth Sci. Chin.* **7**, 255 (2019).

71. Jara-Muñoz, J., Melnick, D. & Strecker, M. R. TerraceM: A MATLAB® tool to analyze marine and lacustrine terraces using high-resolution topography. *Geosphere* **12**, 176–195 (2016).

72. de Gelder, G. et al. Corinth terraces re-visited: Improved paleoshoreline determination using Pleiades-DEMs. *Geotectonic Research in* **97**, 12–14 (2015).

73. Lorscheid, T. & Rovere, A. The indicative meaning calculator-quantification of paleo sea-level relationships by using global wave and tide datasets. *Open Geospatial Data, Software and Standards* **4**, 10 (2019).

74. King, G. C. P., Stein, R. S. & Rundle, J. B. The growth of geological structures by repeated earthquakes 1. Conceptual Framework. *J. Geophys. Res.* **93**, 13307–13318 (1988).

75. Zhou, Y., Parsons, B., Elliott, J. R., Barisin, I. & Walker, R. T. Assessing the ability of Pleiades stereo imagery to determine height changes in earthquakes: A case study for the El Mayor-Cucapah epicentral area. *J. Geophys. Res. [Solid Earth]* **120**, 8793–8808 (2015).

76. Bagnardi, M., González, P. J. & Hooper, A. High-resolution digital elevation model from tri-stereo Pleiades-1 satellite imagery for lava flow volume estimates at Fogo Volcano. *Geophys. Res. Lett.* **43**, 6267–6275 (2016).

77. Almeida, L. P. et al. Deriving high spatial-resolution coastal topography from sub-meter satellite stereo imagery. *Remote Sensing* **11**, 590 (2019).

78. Taylor, J. Introduction to Error Analysis, the Study of Uncertainties in Physical Measurements, 2nd Edition. (ui.adsabs.harvard.edu, 1997).

79. Kendall, R. A., Mitrovica, J. X. & Milne, G. A. On post-glacial sea level–II. Numerical formulation and comparative results on spherically symmetric models. *Geophys. J. Int.* **161**, 679–706 (2005).

80. Milne, G. A. & Mitrovica, J. X. Postglacial sea-level change on a rotating Earth: first results from a gravitationally self-consistent sea-level equation. *Geophys. J. Int.* **126**, F13–F20 (1996).

81. Mitrovica, J. X., Wahr, J., Matsuyama, I. & Paulson, A. The rotational stability of an ice-age earth. *Geophys. J. Int.* **161**, 491–506 (2005).

82. Johnston, P. The effect of spatially non-uniform water loads on prediction of sea-level change. *Geophys. J. Int.* **114**, 615–634 (1993).

83. Milne, G. A., Mitrovica, J. X. & Davis, J. L. Near-field hydro-isostasy: the implementation of a revised sea-level equation. *Geophys. J. Int.* **139**, 464–482 (1999).

84. Lambeck, K., Purcell, A., Johnston, P., Nakada, M. & Yokoyama, Y. Water-load definition in the glacio-hydro-isostatic sea-level equation. *Quat. Sci. Rev.* **22**, 309–318 (2003).

85. Peltier, W. R. & Fairbanks, R. G. Global glacial ice volume and Last Glacial Maximum duration from an extended Barbados sea level record. *Quat. Sci. Rev.* **25**, 3322–3337 (2006).

86. Edwards, R. L. et al. A large drop in Atmospheric 14C/12C and reduced melting in the Younger Dryas, documented with 230Th Ages of corals. *Science* **260**, 962–968 (1993).

87. Abers, G. A. & McCaffrey, R. Active arc-continent collision: earthquakes, gravity anomalies, and fault kinematics in the Huon-Finisterre collision zone, Papua New Guinea. *Tectonics* **13**, 227–245 (1994).

88. Railsback, L. B., Gibbard, P. L., Head, M. J., Voarintsoa, N. R. G. & Toucanne, S. An optimized scheme of lettered marine isotope substages for the last 1.0 million years, and the climatostratigraphic nature of isotope stages and substages. *Quat. Sci. Rev.* **111**, 94–106 (2015).

89. Jouzel, J. et al. Orbital and millennial Antarctic climate variability over the past 800,000 years. *Science* **317**, 793–796 (2007).

90. Lüthi, D. et al. High-resolution carbon dioxide concentration record 650,000–800,000 years before present. *Nature* **453**, 379–382 (2008).

91. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. *Astron. Astrophys. Suppl. Ser.* **428**, 261–285 (2004).

## Acknowledgements

We thank Kathrine Maxwell, Fiona Hibbert and one anonymous reviewer, for their excellent comments and suggestions on earlier versions of this manuscript. We thank the Centre Nationale des Etudes Spatiales (CNES) for the Pleiades imagery, and the German Aerospace Center DLR for TanDEM-X data (Project DEM_GEOL0584). We also thank Pascal Lacroix for help with DSM calculation, Laurie Barrier for the stacked swath profile code, Eric Lewin for help with statistics, Thomas Lorscheid for help with the indicative meaning calculator, and Gerhard Krinner for discussions on paleoclimate. G. G. thanks Robin Lacassin for the initial suggestion to re-visit the Huon terrace sequence with Pleiades topography. G.G. acknowledges the CNES for his postdoctoral scholarship that made this study possible, as well as postdoctoral funding from the IRD and the Manajemen Talenta BRIN fellowship program.

## Author contributions

G.d.G. contributed with conceptualization, methodology, formal analysis, investigation, and writing (original draft & editing). L.H. contributed with conceptualization, discussion, writing (reviewing), validation, and general supervision. A.-M.P. contributed with discussion, writing (reviewing), validation, and modeling supervision. D.F.-B. made the stacked swath profiles and contributed with discussion, writing (reviewing), and validation. T.P. calculated GIA-corrections and contributed with discussion, writing (reviewing), and validation. D.C. contributed with discussion, writing (reviewing), and validation. C.A. contributed with discussion, writing (reviewing), and validation. K.P. contributed with conceptualization, discussion, writing (reviewing), and validation.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-022-00583-7.

**Correspondence** and requests for materials should be addressed to Gino de Gelder.

**Peer review information** *Communications Earth & Environment* thanks Kathrine Maxwell, Fiona Hibbert, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Adam Switzer, Joe Aslin. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022