Distribution of HLG and HLS elevations

During visits to the TKUB site in 2010, 2011 and 2013, and a visit to the TBAT site in 2013, we surveyed HLG elevations on multiple living coral microatolls. On each coral, we surveyed multiple HLG points, following Meltzner and Woodroffe15, and we calculated an average HLG elevation for each microatoll. Some of the microatolls appeared to be connected to the open ocean, whereas others were clearly in ponded settings; hence, our distribution of living HLG elevations includes both open-ocean and ponded microatolls. For those corals that appeared in the field to be open-ocean, we checked their locations on high-resolution imagery to ensure that we did not fail to recognize potential moating, which can be difficult to recognize in the field when water levels are higher. In most cases, the imagery confirmed that the surveyed microatolls were seaward of any potential ponds on the reef, but in a few cases, the microatoll’s setting could not be unambiguously determined. We therefore classified each microatoll as either clearly open-ocean, clearly ponded or possibly ponded.

We wish to determine the indicative meaning of coral HLS immediately after a diedown, as this is the most direct measurement of RSL, and it is the parameter most easily measured in fossil microatoll slabs. What we were able to survey, however, is coral HLG, a number of years after the most recent diedown. Based on a slab through an unponded living microatoll at TKUB38, we determined that the most recent diedown occurred in 2005, coincident with the lowest predicted tides of the 18.61-year cycle, and that the microatolls would have grown up 0.06–0.10 m between then and our surveys in 2010–2013. We therefore subtracted 0.06 m from all HLG elevations surveyed in 2010, 0.07 m from all HLG elevations surveyed in 2011 and 0.10 m from all HLG elevations surveyed in 2013. (The upward growth rate tends to be slower in the first few years after a diedown, so the corrections are not exactly proportional to the time elapsed since the most recent diedown as shown simplistically in Fig. 2.) The spread of coral elevations in Fig. 3 represents the distribution of elevations of HLS immediately following a diedown. HLG elevations in subsequent years would be higher than the elevations shown, by an amount dependent upon the coral growth rate and the time since the most recent diedown.

Calculation of tidal datums

To determine the indicative meaning of coral HLS, we must determine the range of coral HLS elevations relative to tidal datums at each site. To calculate tidal datums, we used the Oregon State University regional tidal inversion for the Indian Ocean region39,40. We extracted the harmonic constituents for each site and used them to calculate mean high water and mean low water (MLW) using formulas from the Manual of Harmonic Constant Reductions41. We note that, because the Belitung region is characterized by diurnal tides, mean high water is equivalent to mean higher high water and MLW is equivalent to mean lower low water (MLLW) at each site. We also determined highest astronomical tide and lowest astronomical tide (LAT) for each site by first computing predicted tide levels every hour over an 18.61-year tidal cycle, and then finding the maximum and minimum elevations. The tidal datums are shown in Fig. 3; note the substantially larger tidal range to the northwest.

Distribution of HLS elevations relative to tidal datums

To tie the surveyed coral elevations into the tidal cycle, we deployed a portable tide-gauge apparatus (a pressure sensor water-level datalogger, with an accompanying barometric pressure datalogger to barometrically correct data recorded by the water-level datalogger) for just over 5 days at TKUB on northwestern Belitung and for just over 2 days at TBAT on southeastern Belitung. Water-level readings were recorded every 10 s and then smoothed to 1-min intervals. We surveyed the base of the water-level datalogger relative to the corals, and we also periodically (several times per day) surveyed the actual water elevation in a calm (but unponded) area, to validate the water-level datalogger readings.

We extracted tidal predictions from the regional tidal inversion for the Indian Ocean region39,40. We adjusted these predictions for local sea-level anomalies estimated on a daily basis from satellite altimetry by the AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic data) group42. We then plotted the tide-gauge readings (and surveyed water elevations) against the adjusted tidal predictions, and we uniformly shifted the vertical reference frame of the entire survey (including the water-level readings and the coral elevations) to minimize the misfit between the recorded water levels and the adjusted tidal predictions. This placed all surveyed corals into a vertical reference frame relative to mean sea level. The resulting coral HLS elevations are plotted alongside the various tidal datums for each site on Fig. 3.

Microatolls at the TBAT site

At the TBAT site on southeastern Belitung, we found a population of microatolls, spread over a minimum distance of 200 m, with each microatoll at a similar elevation as one another. Each of these microatolls was characterized by a high central dome that was surrounded by low middle concentric annuli and high outer annuli (for example, Fig. 4). The morphology of these corals requires growth under changing RSL conditions. The central dome of each microatoll grew during a period when RSL was high; RSL then fell rapidly, killing the upper portions of the corals; RSL then stabilized at a lower elevation, forming a series of low concentric annuli ∼0.6 m higher than present-day analogues; RSL then rose ∼0.6 m in less than a century, allowing the coral to grow upward to 1.2 m higher than modern living corals.

The biggest and best preserved of this population of microatolls (TBAT-F01) had a radius of 4 m, although a significant sector of the coral had irregularities in its outermost annuli, and the radius in that portion of the microatoll was slightly shorter. The outermost part of the microatoll had also cracked and broken away from the central dome, a consequence of the precarious morphology that resulted from its growth pattern, though it was not difficult to fit the inner and outer parts back together. Nonetheless, because of the cracking and irregularities in the outermost annuli, we extracted two radial slabs from this microatoll, and we used each slab to redundantly constrain the site’s RSL history.

For each slab, X-rays were processed and mosaicked together following the guidelines of Meltzner and Woodroffe15. In particular, the brightness and contrast of the X-rays were adjusted to emphasize annual banding, but care was taken to avoid introducing artifacts in the final photomosaics, particularly at the boundary between individual X-rays. The banding visible in the photomosaics was then traced, and the RSL history recorded by the coral was interpreted. The annotated slabs (with photomosaics removed for clarity) are shown in Fig. 5. Full-resolution X-ray photomosaics of each slab, and all original unmodified X-rays, are available from the corresponding author.

Radiocarbon analysis at the TBAT site

A total of eight radiocarbon samples were dated from TBAT-F01. The radiocarbon dates were modelled using the OxCal calibration program43. We applied the Marine13 radiocarbon age calibration curve22, assuming the marine reservoir correction ΔR≈+89 years, based on a ΔR value established from an early 20th century sample from southwestern Borneo23,44. Although there is considerable uncertainty in any ΔR value and its extrapolation spatially and to samples from the mid-Holocene, we can establish that, whatever ΔR was at our sites at the time, it did not vary in a statistically significant way over the lifetimes of our mid-Holocene corals. This observation is crucial, because it allows us to ignore uncertainties in ΔR if we are concerned with only the relative age, or the difference in age, between two corals at the same site.

We use the following argument to demonstrate that ΔR at TBAT did not vary over time. Comparing the unmodelled calibrated radiocarbon dates, assuming for now that ΔR≈+89 years (with zero uncertainty about that assumed value) and accounting for the number of annual growth bands separating the various samples, seven of the eight ages agree at 1σ and all agree at 2σ (Supplementary Table 2). This is consistent with the hypothesis that the reported laboratory errors and the calibration curve correctly describe the uncertainty: 68% of data should agree at 1σ, and 95% should agree at 2σ. This agreement precludes significant variation in ΔR over the 250-year lifetime of TBAT-F01; if the marine reservoir correction varied by more than a few decades over that period, we would not expect such consistency among the unmodelled radiocarbon dates.

Microatolls at the TKUB site

At the TKUB site on northwestern Belitung, no single coral recorded the complete RSL history from ∼6,750 to ∼6,550 cal years BP, but we compiled a RSL history for the period 6,800 to 6,440 cal years BP from five individual microatolls that all grew over a 3-km stretch. Slabs from each of these corals, TKUB-F04, TKUB-F05, TKUB-F16, TKUB-F19 and TKUB-F23, are shown in Fig. 7.

Radiocarbon analysis at the TKUB site

At least two radiocarbon samples were dated from each TKUB coral, and all dates are consistent with their counterparts from the same coral colony at 1σ (Supplementary Table 2). Although ΔR at TKUB may differ from ΔR at TBAT, the consistency among the unmodelled TKUB dates precludes significant variation in ΔR over the lifetime of each coral at the TKUB site.

RSL reconstruction at the TKUB site

The coral growth history based on the TKUB slabs, plotted in Fig. 8, can be divided into three discrete floating chronologies that have been independently radiocarbon dated, but those dates preclude a combined gap of more than ∼21 years between the floating chronologies. These three floating chronologies have been merged together in sequence based on the dating results. TKUB-F04 and TKUB-F05 are the oldest and highest microatolls, and they constitute the oldest floating chronology. They grew at a similar elevation as one another and overlapped in time. TKUB-F16 also started growing at about the same time, but it was ∼0.6 m lower than TKUB-F04 and TKUB-F05. For at least 70 years, it grew with no indication of a diedown or of even being close to HLS. Within two decades after TKUB-F04 and TKUB-F05 died entirely, presumably from sea-level fall, TKUB-F16 recorded its first diedown; it recorded a second diedown, lower than the first, 18 years later. TKUB-F19 began growing at about this time and recorded its first diedown when TKUB-F16 recorded its second diedown. TKUB-F16 and TKUB-F19 continued to grow and to track RSL for nearly a century, forming the middle floating chronology. TKUB-F23, the youngest of the five corals, forms the third floating chronology. Its elevation and growth history suggest that, within two decades after the death of TKUB-F19, RSL rose rapidly, up to a peak only ∼0.2 m lower than the earlier peak recorded by TKUB-F04 and TKUB-F05. More than a century after TKUB-F23 began growing, RSL fell rapidly over less than a decade or two, then gradually rose again over the following ∼30 years.

Types of observations from microatoll slabs

We distinguish four types of observations from a coral slab: uneroded HLS elevations immediately following a diedown; uneroded HLG elevations immediately before a diedown; uneroded HLG elevations in years during which no diedown occurred, when the coral was in unrestricted upgrowth mode; and eroded HLG elevations (the highest level of preserved coral growth) for which it is unknown whether a diedown occurred. The first data type (HLS) is the most direct measurement of RSL, but it tracks only the most extreme low tides and may be biased by an unusual climate or weather event that results in a short duration lowering of sea level. The other data types (HLG) are all technically minimum bounds on low water level, because their elevations are controlled by the coral growth rate and not by RSL. The second data type (HLG just before a diedown) is considered to be a closer approximation to RSL than the third and fourth data types, but such data points are rarely preserved14.

Vertical uncertainty

We distinguish two types of vertical uncertainty in our study. The first is aleatoric and quantifiable: random errors that affect the elevation of one part of a curve relative to another part of the same curve. This accounts for the natural distribution of HLS elevations in any population of corals, including the possible effects of unrecognized ponding. Ponding is a phenomenon whereby some corals can survive at higher elevations than they could otherwise, in elevated enclosed pools that do not drain fully at low tide15,45,46. Ponding is not always easy to recognize, as the effect can be gradual: one pool may raise the water level at extreme low tide by only a few centimetres over the level in an adjacent pool immediately seaward. Nonetheless, the cumulative effect of multiple subtle ponds at progressively higher elevations tends to exceed 0.1 m only on the wider and more physiographically complex reefs46.

To estimate a formal uncertainty about the elevation of any one RSL proxy data point, we surveyed a distribution of HLG elevations on living corals (including some that were clearly ponded) at each site. We augmented this data set with the elevation differences between coeval diedowns seen in slabs from two different living corals at the TKUB site38. The s.d. of differences in elevation of coeval HLG or HLS at each of our Belitung sites is 0.090 m. This is consistent with observations in Australia, but slightly larger than estimates from off the west coast of Sumatra15. The wider distribution of coral HLS on Belitung than off the west coast of Sumatra may occur because of the wider reefs on Belitung, and/or because the tidal range there is larger.

Because ponding is a concern in sea-level studies using coral microatolls, we specifically address whether our results might be biased by ponding in ways that we have not yet considered. At the TKUB site, because the RSL curve was constructed from five separate corals, it is possible that some of the higher and more landward corals (TKUB-F04, TKUB-F05 and/or TKUB-F23) were ponded by significant amounts, that is, by >0.1 m. However, the amplitude of the mid-Holocene oscillations is twice the range of HLS observed among living microatolls on the modern reef, even considering the highest ponded corals (Fig. 3). At the TBAT site, ponding is less likely to explain the observed oscillations, as the oscillations are entirely recorded on individual microatolls. Finally, the two sites are located 80 km apart, on opposite sides of Belitung Island (Fig. 1). This separation is sufficient that it would require a remarkable coincidence to explain the similar changes at the two sites if those changes were caused primarily by localized ponding at each site.

The second type of vertical uncertainty is epistemic and affects the elevation of the entire RSL curve as a whole. These systematic vertical errors are not shown on any figures in this paper, but include uncertainty in the change in tidal range at each site; uncertainties in tectonic effects or compaction at each site; and uncertainty in the HLS elevation of living corals at each site, which is used as the reference elevation for past RSL15,47. These errors are difficult to quantify, but they are likely small. Tide modelling (Supplementary Fig. 2) and tectonic modelling (Supplementary Figs 3–4) suggest both of those effects are on the scale of centimetres, and neither compaction of the thin sediments underlying the fossil corals nor ponding of the living microatolls is likely to bias the RSL curve at a site by more than ∼0.1 m.

Sea-level statistical model

To analyse the RSL proxy data, we constructed an empirical hierarchical statistical model25, separated into three levels: a data level, which models the recording of RSL by proxies; a process level, which models RSL at the different sites; and a hyperparameter level, which characterizes key attributes of the first two levels.

At the data level, RSL index points (HLS elevations following diedowns) from Belitung are preserved typically once or twice per 18.61-year nodal tidal cycle, whereas minimum limiting data (HLG elevations, or minimum bounds on low water level) are resolved each year. We use all of the index points, as they are indicative of sea level. The selection of limiting data is more complicated, however, as our model treats limiting data as faithful sea-level indicators, yet in reality some limiting data are severe underestimates of sea level. Specifically, any limiting data from before a microatoll’s initial diedown represent coral growth up to that initial HLS, and these data may be decimetres (or even metres) below HLS14. Even after a coral’s initial diedown, due to patterns of microatoll growth over the 18.61-year tidal cycle (Fig. 2), some limiting data from our sites are expected to be as much as 0.20 m lower than the theoretical HLS; in these cases, the highest limiting data point within each 18.61-year cycle should be a reasonable approximation of theoretical HLS for that year, and therefore a useful proxy for RSL. In principle, erosion should also be considered at the data level, but because we selected slabs that were well preserved, erosion was negligible (∼0.05 m or less) and can be ignored over most of the time series in our study. An exception to this is the later RSL peak at both sites, ca. 6,600–6,550 years BP, where no diedowns are preserved and erosion may locally exceed 0.15 m. Because of this limitation, our model may underestimate the elevation of the second RSL peak, and the amplitude of the fluctuations we infer in our study should be considered a conservative minimum estimate.

Our preferred strategy for modelling limiting data from the Belitung sites is, therefore, to subsample the limiting data by selecting only the highest limiting point in each 18.61-year bin (Supplementary Fig. 6); nonetheless, we also consider an alternative strategy, in which we use the highest limiting point available for each year (the only point available in most years), excluding only the early part of TKUB-F16, before the coral had grown up to HLS (Supplementary Fig. 7). The preferred strategy is an attempt to use only data that reliably approximate a given year’s theoretical HLS; the alternative strategy is an attempt to use as much of the limiting data as is possibly justifiable.

We model noisy proxy observations (y i ) of RSL elevation as

where i indexes data points and j indexes sites, and the function f j (t) is RSL at site j and time t. Each observation belongs to one of four floating chronologies (the entire record at TBAT, plus three discrete floating chronologies at TKUB), indexed by k ε [0, 3]; each floating chronology is associated with an age shift Δ k . The sea-level observation errors, ɛ i , are treated as uncorrelated and normally distributed, with a s.d. of 0.09 m determined as discussed in the text and earlier in Methods.

Coral ages are constrained by radiocarbon dating methods. Because we can assume that the marine radiocarbon reservoir correction, ΔR, is fixed over time at each site, the relative age uncertainties between the three floating chronologies at the TKUB site are determined by the radiocarbon ages presented in Supplementary Table 2; these inter-slab age uncertainties result in the possibility that one, two or all three of the TKUB floating chronologies are as much as 21 years older. In addition, uncertainty in ΔR at each site allows for an inter-site relative age shift between the overall time series at the TKUB site and that at the TBAT site of up to approximately ±120 years (the ±85-year uncertainty from each site added together in quadrature). Because the modelling depends only upon relative ages and not upon absolute ages, and because the inter-site relative age uncertainty is so much larger than the intra-site relative age uncertainties, we need only three age-shift parameters, {Δ 0 , Δ 1 , Δ 2 }, and we can define them in a way that is more intuitive than elicited by the formula above (we can fix Δ 3 at 0 year). For convenience, we hold the time series at TBAT fixed to that determined assuming ΔR=+89 years, as discussed in the text. Δ 0 is the overall age shift of the TKUB record relative to the TBAT record, and we allow –120 years≤Δ 0 ≤+120 years. Δ 1 and Δ 2 are the age shifts of the oldest and youngest floating chronologies at the TKUB site relative to the central floating chronology at the site, such that the sum of Δ 1 and Δ 2 is a maximum of 21 years (and a minimum of 0 year), where Δ 1 and Δ 2 are shifts in opposite directions, Δ 1 making the oldest slabs older and Δ 2 making the youngest slab younger. Age uncertainties within individual floating chronologies are not incorporated into the model, as the law of superposition prohibits swapping the order of data derived from successive annual bands, effectively rendering the relative age uncertainty to be negligible.

At the process level, f j (t) is specified as the sum of a common (shared) regional sea-level signal g(t), a periodic signal representing the 18.61-year nodal tidal cycle p j (t), a site-specific offset c j , and high-frequency variability w j (t):

The prior distribution of the shared signal, g(t), is a mean-zero Gaussian process48 (GP) characterized by hyperparameters that comprise an amplitude σ g and a timescale of variability τ,

where ρ is the Matérn correlation function with smoothness parameter 3/ 2 and scale τ. The use of a smoothness parameter of 3/ 2 ensures that the first derivative of the process will be defined everywhere, but allows for abrupt changes in rate.

The prior distribution of the periodic signal representing coral growth over the nodal tidal cycle, p j (t), is a mean-zero GP characterized by hyperparameters that comprise an amplitude σ p , a smoothness parameter ν p and a fixed period corresponding to the nodal tidal period, 18.61 years49:

where t and t′ are defined in years. The hyperparameters of this periodic component48 are tuned for each site to simulations of coral growth under present-day nodal tidal cycles at the site. We assumed a coral growth rate r that is normally distributed with a mean of 12 mm per year and s.d. of 2 mm per year and a periodic cycle with tidal amplitudes (σ p ) of 0.186 and 0.089 m at TKUB and TBAT, respectively (Fig. 2). For tuning these hyperparameters, simulated RSL is given by:

The simulated coral height at any given time, CH(t), is equal to the minimum of RSL(t) and the potential growth of the coral according to the randomized growth rate, based on the coral height in the previous year, CH(t–1)+r:

We generate five random, 100-year-long time series at each site and fit these synthetic coral height data to a mean-zero GP, equivalent to the periodic component of the process model plus white noise. We use these maximum-likelihood parameters from this exercise as the amplitude and smoothness hyperparameters in p j (t) of the original process level, above.

The prior distribution of the constant site-specific offset, c j , is normal with mean zero and variance . We restrict this site offset to being constant because we do not expect any physical processes to give rise to significant centennial-scale or sub-centennial-scale variations in RSL between the two sites. The two sites should be exposed to essentially indistinguishable dynamic sea-level changes, and any tectonic deformation at these sites should be small and similar at the two sites (discussed later in Methods). While GIA and changes in tidal range do vary spatially, any changes due to these processes should be small enough on a centennial scale that they are well within any noise.

The prior distribution of the high-frequency variability in RSL, w j (t), is modelled as white noise, with a normal distribution with mean zero and variance , and its additional homoscedastic (equal variability) noise.

We employ an empirical Bayesian analysis method, in which the age-shift parameters {Δ 0 , Δ 1 , Δ 2 } and the hyperparameters {σ g , τ, σ c , σ w } are point estimates calibrated based on the data to maximize the likelihood of the model. The hyperparameters {σ p , ν p } are optimized as described above, based on the present-day tidal cycles and coral growth models at TKUB and TBAT, and are held constant during the optimization of the other hyperparameters. The key output of the model is an estimate of the posterior probability distribution of the RSL field, f j (t), conditional on the tuned hyperparameters (Supplementary Figs 6–9; Supplementary Table 3).

In the end, the model based on our preferred strategy does a reasonable job of separating the non-linear and periodic signals (Supplementary Fig. 6), and the rates of RSL change it estimates should reflect secular trends, minimally biased by vagaries of coral growth variability over the 18.61-year tidal cycle. The alternative model, in contrast, does a poor job of separating out the periodic term, and it forces more high-frequency variability into the non-linear signal, likely overestimating short-term rates of sea-level change. Although we suspect that the high-frequency variability (period ∼30 years) seen only in the alternative model (Supplementary Fig. 7c) is an artifact of that model trying to fit limiting data that severely underestimate theoretical HLS, the fact that both strategies yield fluctuations at a 200-year timescale with peak-to-trough amplitudes of 0.5–0.7 m and similar timing suggests that these model results are robust.

Reinterpretation of published data from southern China

Yu et al.17 surveyed, sampled and dated a suite of coral microatolls from a site on the Leizhou Peninsula, along the southern coast of China. Unlike in our study, where we collected and analysed full radial slabs of each microatoll, they presented primarily point data from the upper surfaces of microatoll annuli. In total, they published 13 dated samples, each of which was tied to the elevation from which it was collected. They also provided photos and cross-sectional sketches of each microatoll, so although those authors focused only on the upper surfaces, they provided enough information to estimate the timing and elevations of the more prominent diedowns.

We reinterpreted the RSL curve of Yu et al.17 (Supplementary Fig. 10) by estimating the timing and elevations of those more prominent diedowns. The reported ages were based on U-Th techniques (typically with small errors) and were all in the expected sequence (ages from the outer annuli of each microatoll were sequentially younger than ages from the inner annuli), so it was straightforward to estimate the timing of each diedown, and to correlate diedowns from one coral to another. Numerous points in each photograph were marked with surveyed elevations, providing a sense of scale, so we were able to estimate the elevations of those diedowns with sufficiently conservative vertical errors.

Last, we correlated coeval annuli from one coral to another based on their reported ages. Again, this was straightforward, as the microatolls provide a consistent, reproducible RSL history, with the same number of prominent diedowns on the various microatolls between any two dates. The advantage of correlating the annuli manifests when considering the handful of U-Th ages in their study that had sizable errors. In the few cases where the chronological errors were so large that the sample age overlapped with sample ages from adjacent annuli, our effort to group the age–elevation data based on the coral morphologies allowed us to minimize the ambiguity of whether a particular sample belonged on one downward swing of the RSL curve or on the subsequent upward swing (Supplementary Fig. 10).

The most notable difference between the Yu et al.17 interpretation and our reinterpretation of the Leizhou Peninsula data is that, in our analysis, the amplitudes of the RSL fluctuations have increased. To give credit to the original authors, Yu et al.17 acknowledged that their RSL curve ‘only represents minimum cycles of fluctuation, because the low-lying gouges were not dated’ and that ‘the amplitudes of sea-level fluctuations should also be treated as representing minimum values’. Yu et al.17 had observed fluctuations based simply upon the differential heights of the sequential annuli (although they missed the third RSL drop, which is required by the third diedown on their microatoll-3), and adding the HLS (diedown) elevations to their RSL curve causes the troughs of the curve to drop lower. The lowering of the RSL troughs is robust even allowing for 0.1–0.2 m coral growth over the 18.61-year nodal tidal cycle, as we have done for our Belitung time series.

We conservatively increase the vertical errors reported by Yu et al.17 to match those determined for the Belitung sites, and we arbitrarily double the vertical errors for the diedowns, since they are estimated from photos. Furthermore, unlike the Belitung time series where data points were annually resolved and relative age uncertainties were negligible, here there are only several limiting data points per century, each has errors spanning two 18.61-year cycles or more, and only the most prominent diedowns (one per century) have been identified. Hence, it is not practical to isolate the 18.61-year periodic component of the time series from the Leizhou Peninsula. Instead, we add additional vertical error (in quadrature) to account for the fact that this signal is not modelled. The reinterpreted data set, with the increased vertical errors and added diedowns, appears in Supplementary Table 7 and Supplementary Fig. 10. Note that Yu et al.17 had dated two samples (FPO-26 and FPO-30) from the same annulus of the same microatoll at the same elevation, and the samples provided similar ages; we combine them here and derive a weighted average of the ages.

Modification of the sea-level model for southern China data

Because the large chronological errors on a few data points in the Leizhou Peninsula data set allow for the swapping of the ordering of those points in ways that are clearly impossible (based on coral morphologies and the law of superposition), it is desirable to trim the chronological errors where appropriate. We ran the original U-Th ages in Supplementary Table 7 through OxCal43, classifying each grouping of data as an unordered Phase() in OxCal; each of these nine phases was separated from adjacent phases by a Boundary(), thereby imposing the sequential ordering of diedowns that is known from the coral morphologies. The OxCal-trimmed ages, which were used in subsequent analyses, appear in Supplementary Table 7 under the heading ‘modelled ages’.

As noted earlier, the Leizhou Peninsula data comprise only a few data points per century, each with errors spanning two 18.61-year cycles or more, rendering it unjustifiable to model the periodic tidal cycle signal; hence, the Leizhou Peninsula model includes only a non-linear signal and a white noise term. The non-linear signal is based on the data from the higher-resolution Belitung sites: it retains the optimized timescale derived from the Belitung sites, but we scale the amplitude hyperparameters by a factor of two compared with the Belitung sites, to compensate for the fact that with the sparser data set, the model with the unscaled amplitude hyperparameters systematically under-predicts the RSL peaks and over-predicts the RSL troughs. An additional complication is that the Leizhou Peninsula data set contains 12 minimum limiting data points (which are clustered near and constrain the peaks of the RSL curve) but only four index points (which are isolated at and constrain the troughs of the RSL curve). An unweighted model tends to fit the limiting data at the expense of fitting the index points; this is effectively a sampling bias. To compensate, we triple the weight on each index point, so that the index points as a whole (the data constraining the troughs) and the limiting data as a whole (the data constraining the peaks) are weighted equally. This has the desired effect that the model fits the peaks and the troughs of the RSL curve more equitably. The model is plotted against data with their original U-Th dates in Supplementary Fig. 10 (this figure also shows the effect of scaling the amplitude hyperparameters), and it is plotted against data with model-refined ages in Fig. 9.

Sea-level model cross-validation

Cross-validation is used to compare the performance of different predictive modelling procedures. For the preferred model (Supplementary Fig. 6), we performed an exhaustive (64 runs, one for each training point) Leave-One-Out Cross-Validation of the model. Since the model is tuned to envelop 95% of the data, we expect the point that is left out of the optimization of the model to be included ∼95% of the time. Supplementary Table 8 shows the number and percentage of data points that were within the 95% interval of our model’s posterior predictive distribution. The model achieved 92.2% inclusion within the prediction interval. In addition, Supplementary Table 8 shows the mean, median and median absolute value of all of the differences (or residuals) between predicted RSL and sea-level height of the data point. For the mean and median, over-predictions and under-predictions tend to cancel one another out, so values near zero suggest that the differences are randomly distributed. For a model that treats each training point as a sea-level index point, we expect such behaviour. The median absolute error is the median of the absolute value of each difference, so values near zero suggest better predictive power of the model.

Differential GIA across Belitung

Modelling GIA is a complex problem and has been the focus of much research in recent decades. There remain large uncertainties both in spatio-temporal details of the ice melting history50 and in the most appropriate rheological model (lithosphere thickness, and upper and lower mantle viscosities) to use. Significant differences persist in the values of these parameters assumed by different global GIA models51,52,53,54. However, a recently developed GIA model for the Southeast Asia region consistently shows that RSL over the past 7 kyr should be higher at site TKUB on northwestern Belitung than at site TBAT on southeastern Belitung, regardless of the choice of earth model and ice model28.

At far-field sites, following ice melting and the inundation of broad, shallow continental shelves, the RSL signal is driven primarily by two GIA processes: equatorial ocean syphoning and continental levering27. Equatorial ocean syphoning results in far-field RSL fall, due to the migration of water from the far field to the near field to fill the regions vacated by the collapsing forebulge. Continental levering from increased ocean load along continental margins induces uplift at inland regions and subsidence within the ocean basin, generating large sea-level gradients perpendicular to the coast (as shown on Supplementary Fig. 1a). It is the interplay between these two processes and the spatially complex signal resulting from the two larger nearby landmasses of Sumatra (to the west) and Borneo (to the east) that drive the difference in the RSL signal between TKUB and TBAT.

In a bathymetrically simple region, such as off the south coast of China, features formed at sea level (such as a wave-cut notch, abrasion platform, or coral reef) during the mid-Holocene on a small island far offshore would now be below sea level, even if sea level itself has not risen since the mid-Holocene. In contrast, similar sea-level markers formed inland (such as in an embayment) would now be higher than when they formed. In the narrow region between Sumatra and Borneo, the GIA signals from these two larger islands interact; this interaction drives a differential uplift signal across Belitung Island (for example, Supplementary Fig. 1a), with TKUB subjected to a larger ocean-load driven RSL fall than TBAT.

Supplementary Fig. 1b shows a χ2 analysis of various rheological models used to predict RSL, from Bradley et al.28. Several potential models are listed in Supplementary Table 5. The ‘96p28’ model, which yields the best fit overall to Holocene data from the Malay–Thai Peninsula, predicts that the RSL at 7 kyr should have been 0.23 m greater at TKUB than at TBAT (Supplementary Table 5 and Supplementary Fig. 1a–d). The ‘9611’ model, which yields only a slightly poorer fit to the Malay–Thai data but produces a marginally better fit to Holocene data from China (but which is still outside the 95% confidence limit for the preferred earth model for China), predicts that the RSL at 7 kyr should have been 0.40 m greater at TKUB than at TBAT (Supplementary Table 5; Supplementary Fig. 1b–d). The ‘96p510’ model would be considered a global-average earth model, and although it falls outside the 95% confidence limit for the preferred earth model for the Malay–Thai Peninsula, it would predict that the RSL at 7 kyr should have been 0.39 m greater at TKUB than at TBAT (Supplementary Table 5; Supplementary Fig. 1b). We note that the earth models that fit the China data well and those that fit the Malay–Thailand data well constitute two generally distinct populations of models. This might be expected, given the significantly different tectonic regimes across the two regions. Although details remain to be resolved and efforts to do so are an active research area, a consistent conclusion from the range of plausible models considered by Bradley et al.28 is that the RSL at 7 kyr was decimetres higher at TKUB than at TBAT. This is supported further by the model of Peltier51, which incorporates the ICE-5G (VM2) ice model–earth model combination and predicts that RSL at 7 kyr was 0.38 m higher at TKUB than at TBAT.

We therefore conclude that much of the 0.5–0.7 m discrepancy between the absolute elevations observed for mid-Holocene RSL at TKUB on northwestern Belitung and TBAT on southeastern Belitung can be explained and is predicted by GIA.

Tidal range across Belitung and its change over time

Belitung Island sits at an exceptional location on a map of tides. As shown on Fig. 3, the tidal range is substantially larger at site TKUB than only 80 km to the southeast at site TBAT. Supplementary Fig. 2a, which shows the spatial distribution of the elevation of MLLW relative to mean sea level, reveals a tight gradient in tidal amplitudes across the island. Previous modelling studies55,56,57 have demonstrated that this gradient is primarily due to spatial differences in the amplitudes of diurnal tidal constituents (K 1 and so on). The variability of the K 1 constituent is a shelf-resonance; the length and width of the basin produce a natural period of oscillation that is closely aligned with the period of the K 1 tide. We wondered whether the resonance pattern or the tidal range at either site might have been different during the mid-Holocene, when local RSL was 1–2 m higher. This is important because HLS tracks lower water levels, and any change over time in the tidal range could bias our reconstructions of RSL during the mid-Holocene.

Answering this question required the application of a dynamical tidal model rather than assimilative model such as TPXO7.2 (ref. 58). Therefore, the two-dimensional depth-integrated version of ADCIRC59 was applied to the region. A large area (Supplementary Fig. 2d) was modelled in order to place the model open boundaries in regions of deep water (Indian Ocean and Pacific Ocean), where small depth changes would not be expected to affect the tidal constituents. ADCIRC uses an unstructured triangular mesh, and mesh resolution was adjusted to place the highest resolution (kilometre scale) in the Java Sea. The final mesh had ∼1 million elements. Mesh bathymetry was drawn from a blend of SRTM30 data60 and global ASTER data. Initial attempts to use ETOPO1 bathymetry61 revealed strongly non-physical artifacts in our study area.

The open boundaries of the model domain were forced with tidal constituents drawn from the FES2004 tidal model62. For the model runs representing higher RSL, the FES2004 constituents were again used. The inherent assumption is that small changes in water depth in the deep ocean locations of the open boundaries will not affect the tides. This assumption is supported by Hill et al.63, who found that large changes in water depth could markedly affect regional or shelf tides (for example, in the western North Atlantic Ocean) but that metre-scale changes had little impact in deep basins. To compute tidal datums, a model run of 90 days was conducted, with the first 30 days as a ramp-up period. Harmonic analysis was performed on the remaining 60 days, yielding amplitudes and phases for each tidal constituent within the domain. These amplitudes and phases were subsequently converted to tidal datums (mean higher high water, MLW and so on) using the Harmonic Constant Datum method64.

Supplementary Fig. 2a–b shows the results for MLLW for the baseline condition and the condition where the water depth has been uniformly increased by 2 m. As discussed earlier, the geometry of the Java Sea is resonant with the K 1 tide and the northwest–southeast gradients in MLLW are essentially a proxy for the spatial variability in the K 1 tide (and, to a lesser extent, in the other diurnal constituents).

Supplementary Fig. 2c shows the change over time of the MLLW elevation. For both sites, the model predicts that MLLW would have been several centimetres lower under the conditions of ∼6.5 kyr, but the change would have been larger at TBAT. For LAT, the changes over time would have been approximately twice those for MLLW: at ∼6.5 kyr, LAT would have been 0.00–0.05 m lower at TKUB on northwestern Belitung and 0.05–0.10 m lower at TBAT on southeastern Belitung.

Two types of modelling artifacts show up in the contour maps (Supplementary Fig. 2a–c) and should be ignored. The occasional crescent ‘artifacts’ observed in the contour maps are due to the region being at a boundary between predominantly semi-diurnal and predominantly diurnal basins. This boundary is quantified by the amplitude ratio R64. The Harmonic Constant Datum method uses slightly different techniques for diurnal and semi-diurnal regions and these slight differences are responsible for the crescent features in Supplementary Fig. 2a–c. In a basin that is strongly semi-diurnal, or strongly diurnal, these artifacts would be absent. In addition, the baseline scenario (Supplementary Fig. 2a) has elements that ‘dry out’ at extremely low tide. The tidal signals at these nodes have truncated troughs, which produce unreliable estimates of amplitudes and phases from the harmonic analysis. These results then propagate into the calculation of tidal datums, and they propagate further into the map of the change over time; the strongly blue bits of Supplementary Fig. 2c are likely an artifact of MLLW being incorrect in those areas in Supplementary Fig. 2a.

These modelling results, taken in consideration of the observed fossil coral elevations, imply that, during the mid-Holocene, both mean sea level and LAT were higher than today, but there was a greater separation between mean sea level and LAT. If the corals directly indicate that LAT at ∼6.5 kyr was at the elevations plotted in Figs 6 and 8 (compared with LAT today), then mean sea level at ∼6.5 kyr would have been 0.00–0.05 m higher at TKUB and 0.05–0.10 m higher at TBAT (compared with mean sea level today) than shown on those plots. These differences are small compared with other errors and can be ignored for most purposes, but we note that a small fraction (perhaps 0.05–0.10 m) of the 0.5–0.7 m discrepancy between the apparent elevations observed for mid-Holocene RSL at TKUB on northwestern Belitung and TBAT on southeastern Belitung can be explained and is predicted by changes in the tidal range over time.

Inferred tectonic stability of Belitung from GPS data

Although Belitung Island is considered tectonically stable, few data exist with which to test any hypothesis of regional tectonic stability. Simons et al.16 published and analysed GPS data from >100 sites across the Southeast Asia region, spanning 10 years, from 1994 to 2004. They defined a relatively undeforming Sundaland block and characterized its boundaries (Fig. 1). Their results suggest that Belitung Island is tectonically stable, at least over the period of their study. While some studies suggest that geodetic deformation can vary on multi-decadal timescales and that 10 years of data are far from sufficient to understand deformation on timescales of centuries to millennia, at least above the seismogenic zone of a subduction megathrust20, there are no mapped active tectonic faults near Belitung, and no evidence has yet been recognized that suggests the Belitung vicinity has been tectonically active over the Holocene.

Potential viscoelastic response to megathrust rupture

The suggestions of tectonic stability in recent decades notwithstanding, Belitung is only 700 km from the Sunda–Java trench, and places such as Phuket, Thailand, at a similar distance from the 2004 rupture, have experienced substantial vertical deformation as part of a viscoelastic response to the 2004 earthquake65. We therefore modelled the potential effects on Belitung of two scenario ruptures along the Sunda megathrust.

We developed models in VISCO1D, using two end-member rheologies previously determined for the Sunda megathrust, to predict the viscoelastic response following a hypothetical rupture along the portion of the Sunda–Java megathrust closest to Belitung Island. The first end-member rheology is that of Pollitz et al.66,67; the second end-member rheology is that of Panet et al.68 Both were constrained by postseismic observations following the 2004 and 2005 Sunda megathrust ruptures, but using different data sets. Parameters of each rheological model are given in Supplementary Table 6.

For each end-member rheology, we modelled the response to two hypothetical ruptures. In each scenario, we assumed the fault ruptures the megathrust up to the surface. We placed the hypothetical ruptures along portions of the Sunda–Java megathrust stretching from southern Sumatra to western Java, which would maximize deformation at Belitung. Although these specific ruptures are neither known nor expected from historical or geological information, the seismogenic potential of this section of the megathrust is poorly understood, and we wished to consider ‘worst-case’ plausible scenarios. The first scenario earthquake for each rheology has M W 8.9, with rupture dimensions of 518 by 175 km and uniform slip of 11.0 m (Supplementary Fig. 3a,e). A rupture of this size in this location may already be pushing or exceeding the limits of what is possible along this section of the megathrust. The second scenario earthquake for each rheology has M W 9.2, such as in 2004, but is more compact, with rupture dimensions of 748 by 236 km and uniform slip of 14.9 m (Supplementary Fig. 3b,f). Rupture dimensions for the chosen magnitudes are based on the scaling relations of Blaser et al.69. The rake angle in both scenarios is 90°, that is, pure thrust motion.

For the rheology of Pollitz et al.66,67, the predicted gradual viscoelastic response at Belitung to a large megathrust rupture centred about the Sunda Strait is trenchward and downward (Supplementary Fig. 3a–d). For the M W 8.9 earthquake, the cumulative vertical displacement at Belitung after 50 years is ∼0.10 m, which is substantially smaller than the amplitude of the oscillations in RSL recorded by the corals. Even for the M W 9.2 rupture, the vertical displacement at Belitung after 50 years is only ∼0.26 m, which is larger but still less than half the amplitude of the observed RSL oscillations.

In contrast, using the rheology of Panet et al.68, the model predicts a gradual viscoelastic response at Belitung that is trenchward and upward (Supplementary Fig. 3e–h). For the M W 8.9 earthquake, the cumulative vertical displacement at Belitung after 50 years is also ∼0.10 m, which again is substantially smaller than the amplitude of the oscillations in RSL recorded by the corals. Even for the M W 9.2 rupture, the predicted vertical displacement at the TKUB and TBAT sites after 50 years is ∼0.28 m or less, still less than half the amplitude of the observed RSL oscillations.

In principle, if the magnitude of viscoelastic deformation were larger (in either direction), it might explain a single oscillation at Belitung. If the postseismic viscoelastic response were downward, then if RSL had reached its peak on the Sunda Shelf and had already begun to fall when, ∼6,750 cal years BP, a large rupture occurred near the Sunda Strait, the response at Belitung over the following 50 years or more would have been one of land-level fall and RSL rise. Once this postseismic viscoelastic response decayed to a negligible rate, RSL fall due to GIA would have resumed, and the RSL highstand at Belitung would have been characterized by a double peak. Alternatively, if the postseismic viscoelastic response were upward, then if the large rupture occurred roughly 150 years before RSL would have otherwise reached its peak due to GIA, then the ongoing sea-level rise would have been overwhelmed by postseismic uplift, and RSL would have begun falling. Once this postseismic viscoelastic response decayed to a negligible rate, RSL rise dominated by the meltwater signal would have resumed, until the RSL highstand, and then after the meltwater production ceased, GIA processes would have led to RSL fall; again, the apparent RSL highstand at Belitung would have been characterized by a double peak.

Nevertheless, we find each of these explanations unsatisfying and unable to explain the totality of the observations for at least four reasons: it would require a coincidence in the timing of the earthquake, either only a few decades after the peak in RSL (in the first scenario) or about 150 years before the peak (in the second scenario); it cannot explain more than two RSL peaks, whereas corals indicate at least three occurred; it cannot explain coeval oscillations 2,600 km away along the south coast of China17 (Figs 1 and 9); and as discussed above, the amplitudes of a predicted viscoelastic response to coseismic rupture simply are not large enough unless we invoke an unrealistically large coseismic rupture to trigger the response or a rheology very different from those that have been published. As a caveat to this last point, we note that the rheologies assumed for the postseismic models (Supplementary Table 6) and the rheologies assumed for the GIA models (Supplementary Table 5) are quite different. The incompatibilities between the postseismic rheologies and the GIA rheologies are common in such studies. This serves, in part, to illustrate how poorly rheology models are constrained. In any case, the four reasons above are sufficient to discount the likelihood that viscoelastic deformation following coseismic rupture along the Sunda–Java megathrust played any role in our observed oscillations on Belitung.

Potential deformation from an unknown upper-plate fault

The 04 June 2015 M W 6.0 Sabah earthquake, 1,400 km northeast of Belitung, occurred on a fault that was previously unrecognized, highlighting the dearth of understanding of the tectonics of this intraplate region. In the hope of considering all possible tectonic explanations for the observed RSL oscillations at Belitung, we also explore hypothetical deformation that might occur if something like the 2015 Sabah rupture occurred on an unknown fault closer to Belitung. As in the previous subsection, we seek merely to answer the question of whether two of the RSL peaks at the two Belitung sites might be explained by coseismic rupture of an upper-plate fault. Specifically, what minimum moment magnitude would be needed for an optimally located and optimally oriented upper-plate rupture to generate ∼0.6 m of uplift or ∼0.6 m of subsidence simultaneously at both Belitung sites?

Using a dip angle similar to the 2015 Sabah rupture (70°) and the scaling relations of Blaser et al.69, we generated a series of synthetic ruptures up to M W 7.6, each with a depth of slip that allows the rupture to propagate to the surface (Supplementary Fig. 4). The smallest rupture that produces uplift or subsidence of ≥0.6 m at a single site has M W 6.8 (Supplementary Fig. 4a); however, the vertical deformation signal is localized over an area with a maximum dimension of ∼30 km, precluding ≥0.6 m of vertical deformation simultaneously at two sites 80 km apart. The smallest rupture that could produce uplift or subsidence of ≥0.6 m simultaneously at two sites 80 km apart has M W 7.6 (Supplementary Fig. 4e); a fault capable of such an earthquake would necessarily be longer than 80 km and have obvious geomorphic expression. In contrast, if any active upper-plate fault exists near Belitung, it must be short enough, with little enough cumulative slip, to have thus far gone unnoticed by geologists. We therefore conclude that rupture of an upper-plate fault is not a viable explanation for the RSL oscillations observed on Belitung between 6,800 and 6,500 cal years BP.

Potential deformation from deeper earthquakes

We discount the likelihood of vertical deformation near Belitung, on the order of ≥0.6 m, from an intermediate-depth or deep-focus earthquake. Although moderate ruptures at depths of ∼600 km or more occurred in 1957 (M W 7.2) and 1963 (M W 7.1) only ∼200 km south of Belitung, and although a larger (M W 7.5) intermediate-depth event occurred in 2007, 350 km to the south70, it is generally observed that intermediate-depth and deep-focus ruptures occur within the subducted slab and tend to produce only centimeter-scale deformation at the surface.

Data availability

Data and modelling codes that have contributed to the reported results are available from the corresponding author upon request.