A sea-level plateau preceding the Marine Isotope Stage 2 minima revealed by Australian sediments

Further understanding of past climate requires a robust estimate of global ice volume fluctuations that in turn rely on accurate global sea-level reconstructions. An advantage of Marine Isotope Stage 2 (MIS 2) is the availability of suitable material for radiocarbon dating to allow comparison of sea-level data with other paleoclimatic proxies. However, the number and accuracy of sea-level records during MIS 2 is currently lacking. Here we present the history of MIS 2 eustatic sea-level change as recorded in the Bonaparte Gulf, northwestern Australia by reconstructing relative sea level and then modeling glacial isostatic adjustment. The isostatically-corrected global sea-level history indicates that sea-level plateaued from 25.9 to 20.4 cal kyr BP (modeled median probability) prior reaching its minimum (19.7 to 19.1 cal kyr BP). Following the plateau, we detect a 10-m global sea-level fall over ~1,000 years and a short duration of the Last Glacial Maximum (global sea-level minimum; 19.7 to 19.1 cal kyr BP). These large changes in ice volume over such a short time indicates that the continental ice sheets never reached their isostatic equilibrium during the Last Glacial Maximum.


Reconstruction of the sedimentary environment as a sea-level indicator. The transect of cores
were recovered from the Bonaparte Gulf, between a water depth of −120 and −67 m during cruise of KH-11-1 in 2011 (Figs 1, 2, and Supplementary Table S1). The range of water depth was selected to provide a continuous RSL record for MIS 2. Following previous work 23 , lagoonal/estuarine, intertidal, and open marine facies are documented in KH11-1 cores based on the abundance of particular benthic foraminifera, such as Ammonia beccarii (Methods and Supplementary Fig. S1), and terrestrial components (as indicated by C/N ratios and Ca/Ti ratios; Supplementary Figs S2 and S3). A. beccarii is a salt-marsh foraminifer useful for reconstructing Holocene sealevel changes 24 . The lagoonal/estuarine and intertidal facies are dominated by sedimentary components consistent with close proximity to the coastline, while the open marine facies consists of primarily marine materials indicated by low C/N ratios (<10) and high Ca/Ti ratios ( Supplementary Figs S2 and S3). The lagoonal/estuarine facies is evident in sediment gravity cores KH11-1-GC10 and KH11-1-GC11 due to a dominantly terrestrial geochemical signature (as indicated by high C/N ratios (>10) and relatively low Ca/Ti ratios; Fig. 2a, Supplementary  Fig. S2), as well as in KH11-1-GC03, GC06, GC07, GC08, and GC09 (Fig. 2b, Supplementary Fig. S2), which additionally contain a high abundance of A. beccarii. In KH11-1-GC14, GC19, and GC13 the intertidal facies is identified by abundant terrestrial organic matter with high C/N ratios (>20) and plant macrofossils in KH11-1-GC14 and GC19. The latter core also contains a prominent sand layer (Fig. 2a, Supplementary Fig. S3). The remaining intervals of each core contain evidence of strong bottom current activity with low C/N ratios (<10) and relatively high Ca/Ti ratios indicating an open marine setting with increased distance from the paleoshoreline. www.nature.com/scientificreports www.nature.com/scientificreports/ Age-depth model reconstruction. In the lagoonal/estuarine facies, marine macrofossils are generally well preserved, associated with clay-rich sediments, and thus appear to be in-situ (Fig. 2c), which is consistent with tidal modeling results indicating low tidal energy during sea-level lowstands (Fig. 3). Conversely, construction of age models within the open marine facies was problematic due to a low number of well-preserved macrofossils and potential reworking of older materials due to high tidal energy during sea-level highstands (Fig. 3). Local radiocarbon reservoir age (ΔR) in the Bonaparte Gulf appears to be negligible because well-preserved marine and terrestrial macrofossils from the same stratigraphic level in KH11-1-GC19 exhibit consistent calendar ages following calibration by Marine13 and IntCal13 25 , respectively (Fig. 2d, Methods, and Supplementary Table S2). paleo-tidal reconstruction. Modeled results for the K1 and M2 components, representing diurnal and semi-diurnal tides respectively, are greatly affected by glacial/interglacial-scale sea-level variability, while the K1, O1, M2, and S2 components are insensitive to such change ( Supplementary Fig. S4). Lowering sea level from 0 m to −120 m results in a much lower tidal range for the K1 and M2 components (Fig. 3). In particular, tidal amplitude approaches zero once sea level drops below −90 m as large swaths of the bay are exposed. O1 and S2 tidal patterns are similar to those of K1 and M2 ( Supplementary Fig. S4).

Discussion
A depth transect of sediment cores is a powerful tool to overcome the limitation of an individual sea-level index point 13 . This approach assumes that two cores with different paleo-water depths at the same age may be explained by a single sea-level curve plus age uncertainty. Uncertainty in paleo-water depth may be constrained by evaluating the consistency between deep-and shallow-water indices ( Supplementary Fig. S5). If a shallower sea-level point cannot be explained, water-depth uncertainty should be revised. We applied this method to our sea-level curve and re-evaluated previous work (Methods and Supplementary Fig. 4b,c). The depth transect indicates that lagoonal/estuarine and intertidal facies occur within +15 m water depth using points from KH11-1-GC03 (point X in Fig. 4b) and GC06 (point Y in Fig. 4b) at 20.5 cal kyr BP. This range is consistent with the present water depth of the estuarine environment in the Bonaparte Gulf ( Supplementary Fig. S6). Here, a water-depth uncertainty of lagoonal/estuarine and intertidal facies of +15 m encompasses contemporaneous sea-level indices.
RSL can be deduced from the ages and depths of significant facies changes in this transect. The age-depth model of K11-1-GC10 indicates that sea level was above −110 m from 25.9 to 23.6 cal kyr BP ( Fig. 4; modeled median age probability), which is consistent with KH11-1-GC11 indicating sea level above −110 m at 25.7 cal kyr BP. Little lithological (massively bedded clay matrix) and geochemical variation suggest only a minor variation of sea level at this time (Figs 2 and 4, Supplementary Fig. S2). While previous work 26 reports a single datum at ca. 22 cal kyr BP with −94 m sea level, this may reflect reworking due to strong bottom currents during sea-level highstands. The age-depth model of KH11-1-GC03 shows that sea level was above −114 m from 21.3 to 20.4 cal kyr BP (Fig. 4b), which is consistent with the lagoonal/estuarine facies observed in KH11-1-GC06. The lagoonal/ intertidal facies in KH11-1-GC03, GC10 and, GC11 indicate a sea-level plateau from 25.9 to 20.4 cal kyr BP, which is shallower than previously reported MIS 2 sea level 27 (Supplementary Fig. S7). RSL above −117 m during 20.8 to 19.7 cal kyr BP is indicated by GC07 and consistent with sea level reconstructed from KH11-1-GC03, GC06, GC08 and GC09 (Fig. 4b). Previous work in the Bonaparte Gulf 13,14 reports a brackish environment at −120 ± 2 m sea level during 19.8 to 19.1 cal kyr BP in core RS176-GC5 (Fig. 1). This agrees with the end of lagoonal/estuarine facies in KH11-1-GC07 (2σ age range). A rapid sea level rise 13,14 indicates that the LGM sea-level minimum terminated at 19.1 cal kyr BP, corresponding to the demise of lagoonal/estuarine facies in KH11-1-GC06. Lagoonal/ www.nature.com/scientificreports www.nature.com/scientificreports/ estuarine facies in KH11-1-GC06 and GC07, combined with the brackish environment in RS176-GC5 13,14 , support the interpretation that the LGM sea-level minimum in the Bonaparte Gulf occurred from 19.7 to 19.1 cal kyr BP (modeled median probability). Sea level then rose to above −98 m by ca. 17 cal kyr BP as indicated by the presence of an intertidal facies in KH11-1-GC13 (Fig. 4a,c). The intertidal facies observed in KH11-1-GC14 and GC19 further indicate a rapid sea-level rise at ca. 14 cal kyr BP, corresponding to Meltwater Pulse 1A interval 28 .
Our constraints on paleo-tidal amplitude indicate that tidal amplitude is reduced below −70 m sea level due to a shift of an amphidromic area (i.e., where tidal amplitudes take minimum values) and the changes in the basin configuration of the Bonaparte Gulf (Fig. 3). The Bonaparte Gulf was largely isolated from the open ocean due to lowered sea level with restricted connections between exposed carbonate platforms 16,17 (Fig. 1). The location of submerged platforms could also be identified in the −70 m model results because they roughly coincide with a band of K1 amplitude minimum (Fig. 3) or an offshore limit of M2 amplitudes lower than 0.5 m (Fig. 3). The semi-enclosed environment during sea-level lowstands below −90 m could enhance deposition due to the decreased tidal influences in the Bonaparte Gulf and lack of bottom current scouring, suggesting that the environment is suitable to reconstruct RSL change during the exposure of carbonate platforms and terraces.
The 500-km width of the Bonaparte Gulf continental shelf causes a sea-level gradient across this shelf to appear due to a local difference in isostatic effect 13 and bathymetry ( Supplementary Fig. S8), which can be calculated using a GIA model. The facies-based observational RSL records described above can then be corrected to sea-level equivalent values at KH11-1-GC06 (Fig. 4) to consider the offset of GIA predictions among sites 29,30 . This global ice volume history model consistently explains observational RSL, indicating that a sea-level plateau occurred from 25.9 to 20.4 cal kyr BP, prior to the LGM (19.7 to 19.1 cal kyr BP; Fig. 4, Supplementary Fig. S8).
Considering the calculated regional isostatic effect, the Bonaparte Gulf RSL can be compared with other RSL records. Observational RSL records during MIS 2 and 3 using uplifted terraces in the Huon Peninsula 10-12 are consistent with predicted RSL of this study's ice model (Fig. 4d), supporting the previous-suggested uplift rate of the Huon Peninsula. A MIS-2 RSL record was previously created through sedimentary environmental reconstruction and radiocarbon dating of near-shore cores on the Sunda Shelf 8,9 (Fig. 4e). Because the radiocarbon dating was performed on organic matter with leachable and insoluble components, it exhibits ca. 1,000 years difference within the same horizon, likely due to the effect of old carbon 31 . The Sunda Shelf record is also consistent with our global ice volume history within age uncertainties derived from the effect of old carbon 31 . Additionally, MIS 2 RSL has been reconstructed through precise U-series dating of Barbados corals (Fig. 4f) 6,7 . However, inconsistent results prior to 19 cal kyr BP suggest differential uplift rates between sites, perhaps caused by faulting 32 . There is also large uncertainty in the growth positions (20-50 m) of particular species (Montastrea annularis, P. asteroids, and Diploria) 7 .
Rapid ice-sheet growth of approximately 10 m ice-volume equivalent sea-level between the sea-level plateau prior to the LGM and the time of minimum RSL is observed in our RSL records and the GIA modeling (Fig. 4a,b). www.nature.com/scientificreports www.nature.com/scientificreports/ This would primarily be ice in the Northern hemisphere because the maximum volume added to the Antarctic ice sheet during the LGM is reported to be less than 10 m 33 . The extremely short duration of the global ice volume maximum, revealed by our results, has implications for the usage of LGM boundary conditions for modeling MIS 2 climate since continental ice sheets likely never reached isostatic equilibrium. This could lead to overestimating the maximum volume of each continental ice sheet, which in total may be closer to that inferred from far-field sea-level records 34 . However, additional MIS 2 records are needed to more fully explore these implications. In particular, near-field sites with highly-resolved age dating are required to understand the discrepancy between  13,14 . However, this core was not reevaluated in ref. 14 . Other sea-level data from corals 6 and sediments 8 show the sea-level position above −80 m, suggesting that the paleo-water depth of this brackish facies is underestimated. (d) RSL reconstruction from the Huon Peninsula [10][11][12] . (e) RSL reconstruction from the Sunda Shelf 8,9 . Gray shades show radiocarbon ages of acid-insoluble and leachable organic matter in the same horizon. Uncertainties of age derive from the effect of old carbon 30 . (f) RSL reconstruction from Barbados 6,7 . www.nature.com/scientificreports www.nature.com/scientificreports/ reconstructed global ice volume and the total estimated from individual ice-sheets by near-field dataset 34 . We conclude that a sea-level plateau occurred from 25.9 to 20.4 cal kyr BP with a subsequent 10 m sea-level fall in ~1,000 years, suggesting that continental ice sheets during MIS 2 were less stable than previously believed 27 ( Supplementary Fig. S7). Future work reconstructing paleoclimate during glacial periods should consider the instability in sea level and ice volume revealed by this research.

Methods
Geochemical analyses were performed on sediments to determine the relative contribution of marine and terrestrial components to detect changes in the distance from the coastline. C/N ratios and δ 13 C were measured at a 4 to 10-cm interval using EA-IRMS (Elemental Analysis-Isotope Ratio Mass Spectrometry; Flash EA 1112 and Delta plus Advantage) at the Center for Advanced Marine Core Research (CAMCR) in Kochi University 15,17 . X-ray fluorescence (XRF) core scanning was conducted at a 1-cm interval using TATSCAN-F2 at the CAMCR 35 .
Approximately 100 radiocarbon measurements were performed on marine macrofossils, benthic foraminifera, and terrestrial macrofossils (Supplementary Table S2). Marine macrofossils and foraminifera were cleaned using an ultrasonic cleaner and then etched by 10 M HCl to remove secondary and/or contaminating carbonate materials, which is often attached to the surface of samples. Graphitization was performed following previously reported methods 36 . The analyses were performed at The University of Tokyo using a Tandem Accelerator at the Micro Analysis Laboratory and a Single Stage AMS at Atmosphere and Ocean Research Institute 37 .
Radiocarbon ages were calibrated to calendar ages using MatCal 38 with Marine13 for marine macrofossils and foraminifera and IntCal13 for terrestrial macrofossils 25 . The local reservoir correction is not known for the Bonaparte Gulf but is likely very small 39,40 . Therefore we did not apply a local reservoir (ΔR) correction. Age modeling for cores was performed using 75 of the radiocarbon dates in a 2000-iteration Monte Carlo simulation adapted from ref. 41 , where one probability-weighted age per iteration was selected from the 95.4% range of each calibrated age probability density function (PDF). Over 20 radiocarbon dates are excluded from age modeling due to the reversed age-depth relation. Sampling was repeated up to 1000 times if an age reversal occurred. The calendar year age at each stratigraphic level is the median age of all non-reversing iterations, and uncertainty is reported as the 2.3 and 97.7 percentiles (95.4% range).
Paleotidal variations were reconstructed using a two-dimensional model 22 covering longitude 115°E-143°E and latitude 8.5°S-22°S with the resolution of 1/4 degrees. The simulation was conducted by lowering the sea level uniformly every 10 m from 0 m to −120 m and neglected the impact of isostatic deformation on the seafloor because the Bonaparte Gulf is tectonically stable and far from the former ice sheet [13][14][15]42 . Tides were forced by specifying tidal elevations along open boundaries compiled from harmonic constants of diurnal (K1 and O1) and semi-diurnal (M2 and S2) constituents of a modern global tidal model (TPXO8-atlas) 43 . The M2 constituent is the largest semi-diurnal tidal component and K1 is the largest diurnal component, indicating that the variations of these two constituents reflect the general tidal evolution in the Bonaparte Gulf. We applied the present-day values at the model grid boundaries because global paleotidal models show that tidal elevations along the offshore boundaries did not significantly fluctuate during MIS 2 21,44 and also because modern ocean tides are well-constrained by satellite altimetry, which is not available for the geological past.
A reconstruction of the paleoecological environment of Bonaparte Gulf from vibracore data was presented by ref. 23 . In this previous study, six environments were identified based on biological assemblage: intertidal, lagoonal, estuarine, strandlines, shelf, and riverine with occurrence of A. beccarii indicating lagoonal and estuarine environments. Studies using sediments from Europe, Australia, and East China Sea also support lagoonal and estuarine habitats for A. beccarii [45][46][47] . We conducted a qualitative analysis of A. beccarii abundance, classifying it as barren, low, and high ( Supplementary Fig. S1).
The previous work 13,42 assigned an uncertainty of +5 m to a marginal marine facies based on an assemblage indicative of normal sea-water salinity 14 . However, data at ca. 16.5 cal kyr BP ( Fig. 4c; point α) cannot be explained by a single sea-level curve deduced from point β (Fig. 4c), suggesting the need for revision. Studies on RSL during the deglaciation agree with a sea-level curve deduced from point β, requiring adjusting depth uncertainty for this facies to +25 m. This adjustment is also supported by an apparent mismatch in the age and magnitude of the post-LGM sea level rise 29 indicated by other sea-level records such as Barbados 6,7 and Sunda Shelf 8,9 .
The GIA model is based on the below sea-level equation 48,49 . RSL at site ϕ and time t ( ( , t) rsl Δζ ϕ ) are expressed as follows: Δζ is ice volume equivalent sea level (ESL) and ( , t) ice iso δζ ϕ and δζ ϕ ( , t) water iso is the glacio-and hydroisostatic contributions, respectively. Variations in δζ ϕ ( , t) ice iso and δζ ϕ ( , t) water iso are controlled by the earth model, which is described by the elastic thickness of the lithosphere (H), upper-mantle viscosity (η um ), and lower mantle viscosity (η lm ). The parameters of the earth model are based on the PREM 50 for the density and elastic constants. The adopted earth parameters in this study are based on the previous work 27, 31 . We set H, η um , and η lm to 70 km, (1 − 9) × 10 20 Pa s, and (0.5 − 0.9, 1.0 − 10) × 10 22 Pa s, respectively.
The ANU model cannot explain the sea-level plateau from 25.9 to 20.4 cal kyr BP within the range of predictions derived by the earth model ( Supplementary Fig. S7c). ESL in this study's ice model is based on the ANU model and multiplied with different constant factors for each time-step. The ratio of each continental ice-sheet volume in this study's ice model is the same as the ANU model.