Introduction

Linking past sea-level changes to climatic variations and ice-sheet configurations is fundamental to predicting future climate change scenarios1 and understanding large-scale geodynamics2. The Quaternary represents a key geological archive, as it is the most recent period recording the full spectrum of glacial-interglacial sea-level changes3. Continuous records of past global mean sea level (GMSL) are usually derived from oxygen isotope ratios (δ18O) in marine sediment cores, as the seawater ratio between the lighter (16O) and heavier (18O) isotopes relies on the amount of ice stored within continental ice sheets. Since temperatures and local hydrological/chemical conditions also affect δ18O ratios, many calibration techniques have been proposed for the derivation of GMSL curves from δ18O-records, which has led to a broad range of GMSL estimates4,5. Commonly, to verify or adjust δ18O-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 sparse6,7, rely on uncertain corrections for GIA or other processes of vertical displacement, and are mostly limited to interglacial periods8 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)9. 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 δ18O-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−1 in the NW to 3.5 mm yr−1 in the SE, have led to the preservation of an exceptionally complete sea-level record10. 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 ka11).

Fig. 1: Uplift rates, dated samples and tectonic context of the Huon coral reef terrace sequence.
figure 1

a Slope map of coral reef terrace sequence overlain with colored elevations, showing a comparison of the uplift rates previously proposed10 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 sequences87 ranging from mb = 5.5 (smallest) to 6.2 (largest). GPS arrows25 (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 compilation7. All samples in blue, and in orange the samples after removing duplicates, samples out of the study area and U/Th open system samples.

Our data provides a foundation for improved constraints on Quaternary sea-level oscillations and ice-sheets, from which we conclude that δ18O-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 oscillations10,12,13,14,15,16. 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 constrained6,17,18. 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. 110), whereas later studies simply used an uplift rate calculated from the nearest MIS 5e terrace15,16,19. 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 profiles20,21 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 strike22 (hereafter termed stacked swaths23). 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 highstands19. 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-RSL24. 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 constraints25 (Fig. 1). Theoretically, microplate rotation of 8° per Ma26 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.

Fig. 2: Best-fitting tilt direction of Huon coral reef terrace sequence.
figure 2

a Field-observations re-drawn from a previous study10, 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 analysis72. 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 X2 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).

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 work11. 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−1 uplift rate should give the RSL elevation during terrace formation. We performed this calculation (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 terraces10,11. Instead of using a nomenclature based on reef structure10, 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 terrace11 and RSL elevation of 7 ± 5 m. Within the considerable debate on the elevation and timing of MIS 5e sea-level at different locations17,27,28,29,30 we find the range of 2 m to 12 m a reasonable compromise between the lower end of early MIS 5e GMSL estimates29,31, if GIA effects are negligible at Huon, and the upper end of early MIS 5e GMSL estimates17,32, if GIA corrections at Huon are on the order of +4 m27. Some RSL lowstand elevations were previously calculated from the Tewai section (Fig. 114,33), which we updated using our new uplift rate estimates.

Fig. 3: Huon relative sea-level estimates.
figure 3

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 estimates11, global mean sea-level estimates of MIS 335, MIS 5a36, and MIS 5c36, and a global mean sea-level curve obtained with a Principal Component Analysis4 with 97.5% confidence interval (grey band). Lowstand estimates (Huon RSL - Tewai) were re-calculated from other studies11,33 (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 estimates32 (MIS 5e, 11c). Marine Isotope Stages (MIS; bold) follow a standardized definition88.

New RSL estimates for the last glacial cycle, shown in Fig. 3c, are generally higher in elevation than previous Huon estimates11. 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 level34,35,36. 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 ages11, 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. 3c4. Also our re-calculated 5b and 5d lowstands are younger and at higher elevation than previously proposed for Huon11,14, more in line with the displayed GMSL curve4 (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 Lisiecki4, 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 GMSL and Huon RSL, and we highlight that the similarity of the Huon MIS 11c RSL highstand to independent GMSL estimates37,38,39 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.32, 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 δ18O-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 δ18O-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 curve4. 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 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 sequences40. 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 reefs41,42. 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 code42,43. We test the impact of our RSL estimates by considering 3 different sea-level curves (Fig. 4a; Supplementary Fig. 7): (1) a representative GMSL-curve4, (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 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.44), 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;19,20,45).

Fig. 4: Huon terrace sequence modeling.
figure 4

a Sea-level curves used for terrace modeling, including a GMSL curve based on Principal Component Analysis4 (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−1, erosional potential of 30 mm yr−1 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 profile45, whereas the dated samples for the Kanzarua section15,45,58 were taken from a location ~2 km S of the stacked swath profile.

Considering the misfits between modeled and observed terrace elevations, curves with adjusted highstands better fit the observed terrace morphology, particularly for interstadial terraces T11, 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-stratigraphy16 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 terrace46 (Supplementary Fig. 8). The second potential limitation is that a terrace surface might form several meters below the peak in sea-level (SL-TerraceDIF), depending on the rate of sea-level changes and on the morphology formed during preceding sea-level oscillations or erosion of a barrier reef42. We quantified these effects (Supplementary Fig. 8) using the range of uplift rates determined for Huon (0.5–3.5 mm yr−1) 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-TerraceDIF 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 δ18O-derived GMSL estimates appear to be cyclic, rather than increasing gradually (Fig. 3d), neither peak-delays nor important values of SL-TerraceDIF can easily explain the highstand discrepancy for interstadial sea-levels.

Discussion

Collectively, our geometric analysis, and reef modeling results suggest that δ18O-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 rates47,48. 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 Pleistocene49, 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 δ18O-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 δ18O records during MIS 335,50,51,52. Similar trends may hold for the older interstadials peaks within MIS 6, 9a, 11a. Previous studies highlighted several problems in the conversion from δ18O to sea-level4,53,54,55. Along the lines of those studies, we envision two possibilities, either (1) δ18O has been incorrectly calibrated to sea-level in previous work, or (2) δ18O 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.55 proposed that δ18O 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 δ18O-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 δ18O with geologic indicators. An alternative explanation for relatively low sea-level estimates in δ18O-derived curves is the possibility that sea-level during the last glacial maximum (LGM), to which most δ18O-derived curves are calibrated, may be generally estimated too low56.

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 δ18O records. Such records are often stacked or averaged, ignoring potentially important differences in timing between oceanic basins55,57, and smoothing out rapid sea-level changes4. Rapid (7 ky periodicity), high-amplitude (10–15 m) oscillations have been proposed for MIS 3 in Huon16,58 in relation to Heinrich events59, 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 360, 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 modeling16,58 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, CO2, 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.461. After each interglacial peak, temperature and, to a lesser extent, CO2, 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 rebound62,63,64, differences in glacial thresholds of the different ice-sheets65,66, and cooling of the globe as a precondition for ice sheet growth67. Whichever of these mechanisms is dominant to the relationship between GMSL and paleotemperature/CO2 (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.67 showed, using a global stack of paired sea surface temperature-planktonic δ18O records, that rates of GMSL-changes are largely accounted for by the combination of N-hemisphere summer insolation, sea-surface temperatures (SST; 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 δ18O.

Fig. 5: Comparison paleo sea-level and paleo climate data.
figure 5

a Comparison of GMSL4, Huon RSL, temperature89 and CO290. Temperature anomaly and CO2 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 2191 and sea-surface temperature variations from a global stack67. 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 work67) based on the summer insolation and global sea-surface temperatures in b, and sea-level curves in a.

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 Stereopipeline68 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 Topotoolbox 2.069 and stacked them orthogonal to their trend using an in-house script. This wide swath22, later termed stacked swaths23, 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.20 and age interpretations of Esat et al.19 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 TerraceM70 (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.71). Assuming a broad range of 60–90° as paleocliff slope, we used a free-cliff analysis72 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 subdivisions10. 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}$$
(1)

where HSA is the shoreline angle elevation, HRSL the relative sea-level elevation with respect to present-day sea-level, HREF 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 }_{{H}_{{SA}}}}^{2}+{{\sigma }_{{H}_{{RSL}}}}^{2}{{+\sigma }_{{H}_{{REF}}}}^{2}+{{\sigma }_{{H}_{{SEIS}}}}^{2}+{{\sigma }_{{H}_{{DSM}}}}^{2}}{{\left({H}_{{SA}}-{H}_{{RSL}}+{H}_{{REF}}\right)}^{2}}\right)+\left(\frac{{{\sigma }_{T}}^{2}}{{T}^{2}}\right)\right)$$
(2)

where σHSA, σHRSL, σHREF, σHSEIS, σHDSM, and σ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 σ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 HRSL and σHRSL we use 7 ± 5 m here, as a compromise between reasonable MIS 5e minimum and maximum estimates of 2 m3 and 12 m27, respectively. For HREF and σHREF we use 0.66 ± 0.33 m, as calculated with the indicative meaning calculator73 (Supplementary Table 2). For σHSEIS we used the results of Ota and Chappell49 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−1 uplift rate. Assuming that long-term uplift is a cumulative effect of several seismic cycles74, and thus the long-term uplift rate is proportional to the seismic cycle deviation, we use σHSEIS = U/2 as a rough estimation. For σ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 m75,76,77. For T and σT we use 127 ± 2 ka for MIS 5e, as was used in Lambeck and Chappell11 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 Taylor78, we introduced the weights w = 1/σU2, with best estimates for constants A and B calculated from:

$$A=\,\frac{\sum w{l}^{2}\,\sum {wU}-\sum {wl}\,\sum {wlU}}{\Delta }$$
(3)
$$B=\,\frac{\sum w\sum {wlU}-\sum {wl}\,\sum {wU}}{\Delta }$$
(4)

with

$$\Delta =\sum w\sum w{l}^{2}-{(\sum {wl})}^{2}$$
(5)

And the uncertainties in A and B as:

$${\sigma }_{A}=\sqrt{\frac{\sum w{l}^{2}}{\Delta }}$$
(6)
$${\sigma }_{B}=\sqrt{\frac{\sum w}{\Delta }}$$
(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−1, 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−1 (Supplementary Fig. 4). We then added HREF of 0.66 m to every shoreline angle estimate to obtain past sea-level elevations SLP, with a standard error calculated from σHSA, σHREF, σHSEIS and σHDSM. Using Formulas 3–7, but replacing U with SLP, 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:

$${{{{{{\rm{T}}}}}}}_{{TER}}={T}_{T16}\cdot \frac{{B}_{{TER}}}{{B}_{T16}}$$
(8)

with standard error

$${{{{{{\rm{\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}}$$
(9)

Where TT16, BT16 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

$${{{{{{\rm{H}}}}}}}_{{RSL}}={A}_{{TER}}+({B}_{{TER}}\cdot {l}_{u=0})$$
(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}}$$
(11)

In which ATER and BTER are constants to fit each terrace. Finally, to re-calculate the lowstand elevations given in previous work11,14,33 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−1 (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.79 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 level80,81, evolving shorelines and the migration of grounded portions of marine-based ice sheets79,82,83,84.

Our numerical predictions required models for Earth’s viscoelastic structure and the history of global ice cover. We adopted the viscosity profile VM285, 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 × 1020 Pa s, and 5 × 1021 Pa s, respectively, and (2) a lithospheric thickness of 48 km, an upper and lower mantle viscosity of 5 × 1020 Pa s, and 15 × 1021 Pa s, respectively. These Earth structure parameters are in agreement with the range estimated in Lambeck and Chappell11. We created a global ice sheet history associated with the GMSL history inferred in Spratt and Lisiecki4. 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-5G85. 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 Lisiecki4 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 Lisiecki4. 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 Lisiecki4. For the reef modeling, we used a potential reef growth rate of 10 mm yr−1, as was used by Chappell16 based on U-series dating of the Late-Glacial reef86. We used a basal slope similar to the average slope along the modelled profiles, and erosional potential of 30 mm yr−1 (see full description of the parameterization in Pastier et al.42).

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-TerraceDIF; Supplementary Fig. 8). As for Fig. 4, we modelled terrace sequences for an uplift rate range of 0.5–3.5 mm yr−1, potential reef growth rate of 10 mm yr−1, erosional potential of 30 mm yr−1 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-TerraceDIF, 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-TerraceDIF could have influenced our Huon RSL estimates (Supplementary Fig. 8e).