Fraser Island (K'gari) and initiation of the Great Barrier Reef linked by Middle Pleistocene sea-level change

The eastern Australia coastline is characterized by impressive coastal landforms and an extensive northward-moving longshore drift system that have been influenced by a stable, long-term tectonic history over the Quaternary period. However, the timing and drivers of the formation of two conspicuous landscape features—Fraser Island (K’gari) and the Great Barrier Reef—remain poorly understood. Here we use optically stimulated luminescence and palaeomagnetic dating to constrain the formation of the extensive dunes that make up Fraser Island, the world’s largest sand island, and adjacent Cooloola Sand Mass in southeastern Queensland. We find that both formed between 1.2 Ma and 0.7 Ma, during a global climate reconfiguration across the Middle Pleistocene transition. They formed as a direct result of increased amplitude of sea-level fluctuations associated with increasing global ice volume that redistributed previously stored sediment across the continental shelf. The development of Fraser Island dramatically reduced sediment supply to the continental shelf north of the island. This facilitated widespread coral reef formation in the southern and central Great Barrier Reef and was a necessary precondition for its development. This major reorganization of the coastal sedimentary system is probably not unique to eastern Australia and should be investigated in other passive-margin coastlines. Disruption of sediment flows along the eastern Australia coast due to the Middle Pleistocene formation of Fraser Island set the stage for Great Barrier Reef initiation, according to optically stimulated luminescence and palaeomagnetic dating of sand dunes.

The eastern Australia coastline is characterized by impressive coastal landforms and an extensive northward-moving longshore drift system that have been influenced by a stable, long-term tectonic history over the Quaternary period. However, the timing and drivers of the formation of two conspicuous landscape features-Fraser Island (K'gari) and the Great Barrier Reef-remain poorly understood. Here we use optically stimulated luminescence and palaeomagnetic dating to constrain the formation of the extensive dunes that make up Fraser Island, the world's largest sand island, and adjacent Cooloola Sand Mass in southeastern Queensland. We find that both formed between 1.2 Ma and 0.7 Ma, during a global climate reconfiguration across the Middle Pleistocene transition. They formed as a direct result of increased amplitude of sea-level fluctuations associated with increasing global ice volume that redistributed previously stored sediment across the continental shelf. The development of Fraser Island dramatically reduced sediment supply to the continental shelf north of the island. This facilitated widespread coral reef formation in the southern and central Great Barrier Reef and was a necessary precondition for its development. This major reorganization of the coastal sedimentary system is probably not unique to eastern Australia and should be investigated in other passive-margin coastlines.
The Great Barrier Reef (GBR) is regarded as one of the most notable global biodiversity hotspots and carbon sinks 1 , yet initiation of the world's largest coral reef is poorly constrained and the mechanism responsible is not known 2 . Earliest geological reef development phases along northeastern Australia have been dated to the early to middle Miocene (~23-12 million years ago (Ma)), but age control from in situ coral remains suggests that the extant GBR was established much later, during the Middle Pleistocene, by ~450 ka (refs. 3,4 ). Further evidence from distal offshore sediment records 5,6 in south-central GBR and reef boreholes 2,7,8 from the northern GBR suggest initiation by Marine Article https://doi.org/10.1038/s41561-022-01062-6 Data Fig. 8 and Extended Data Table 2), which indicate that the dune fields were established before the Matuyama/Brunhes boundary (773 ka) 12,13 . Before their formation, northward sand transportation along the eastern Australia coastline continued unabated, creating conditions too sand rich and turbid for widespread reef growth. Once Fraser Island and the dune fields were emplaced, they acted as a barrier to longshore drift that redirected sediment off the edge of the continental shelf 14 , facilitating an increase in carbonate sedimentation and reef growth north of Fraser Island.
The Cooloola Sand Mass and adjacent Fraser Island (in the process of being renamed to its Aboriginal name, K'gari), which is the world's largest sand island and a UNESCO World Heritage area 15 , form part of the extensive coastal dune sequences of southeastern Queensland that also include Moreton and Stradbroke islands to the south 16,17 .
Isotope Stage (MIS) 15 (621-563 ka). The relatively young age of GBR initiation contrasts with the long-term tectonic quiescence and stable coastal configuration along eastern Australia, which was established with the onset of the Antarctic Circumpolar Current by 25 Ma (refs. 9,10 ). Climatic conditions along the central and southern Queensland coastline have been appropriate for reef growth since at least the Pliocene (5.3 Ma-2.6 Ma) 11 , which makes the geologically recent establishment of the modern GBR a conundrum.
In this Article, we present a chronology for Fraser Island (K'gari) and the adjacent Cooloola Sand Mass, which consist of a thick coastal transgressive dune succession, and provide evidence that their formation was a necessary precondition for initiation of the southern and central GBR. The earliest dune-building phase is dated to ~1.2 Ma (Tables 1 and 2) and contains reversed-polarity sediments (Extended  Table 2). f These samples were collected above the cliff from exposed units in the Carlo Blowout. Originally published in ref. 22 Article https://doi.org/10.1038/s41561-022-01062-6 Fraser Island is the terminus of one of the world's longest down-drift coastal systems 14 (Fig. 1). Detrital zircon studies indicate that the sand originates from the Lachlan and New England orogens in central New South Wales along with the Triassic Hawkesbury Sandstone of the Sydney Basin and that it travels 1,500 km north along the east Australia coast 14,18 . The net sand volume transported is estimated to be 500,000 m 3 yr -1 (ref. 14 ). This northward longshore transport follows contours along the east coast of Fraser Island and ultimately deposits sand off the continental shelf north of Breaksea Spit ( Fig. 1) 14 . Existing chronologies from the Cooloola Sand Mass indicate that sand has been accumulating since at least 730 ka (refs. 19,20 ). In this study, we targeted basal units exposed along the Rainbow Beach and Fraser Island cliffs (Supplementary Section 1.0 and Extended Data Fig. 1) that extend their earliest deposition to 1 Ma. Our chronology indicates that the largest and most extensive dune packages formed during post-glacial sea-level transgressions [20][21][22] Fig. 1). Age control for the dune sands was obtained using optically stimulated luminescence (OSL) dating of quartz sand 24,25 , which provides an age estimate for the last time sediments were exposed to light and is used to determine past dune activity. The OSL ages from the oldest basal dune packages indicate initial deposition at ~0.8-1.2 Ma (Fig. 2, Table 1 and Extended Data Figs. [2][3][4][5]. These exceptionally old ages are possible because of the low radioactivity of the ultra-mature quartz-rich sediment 14,18 and have been replicated by other studies 19,20 and validated through multiple samples within the same units (Table 1). While there is OSL age scatter for the oldest dune units, the reported uncertainties are relatively low (10-12% relative standard error) and within the expected range for OSL dating 26 . The accuracy of the basal OSL ages is also supported by the reversed or transitional-polarity   Stratigraphic and sedimentological observations along the Rainbow Beach cliffs indicate two major dune morphologies. The modern Cooloola Sand Mass is composed of large parabolic dunes represented stratigraphically by relatively thin (<10 m) cross-bedded sand packages in the upper cliff section. These units have laterally variable thickness and represent dune emplacement in a vegetated landscape. By contrast, the lower three units exposed in the cliff section comprise 10-to 30-m-thick tabular sets of cross-bedded dune sands with more planar bedding surfaces (Extended Data Figs. 2-4). These packages are traceable for up to 8 km along the Rainbow Beach cliff section. These Middle Pleistocene dune-forming events are interpreted to represent large-scale transgressive dune-field activation under a high sediment-supply regime. Aeolian units with similar age and morphostratigraphy were identified on Fraser Island (Extended Data Fig. 5), which indicates that both dune fields were part of the same Middle Pleistocene coastal dune system 16,22 .
The 700 ka to >1 Ma ages and reversed polarity in basal aeolian units suggest that the Cooloola and Fraser Island dune fields were established during the MPT. Earth's climate system changed fundamentally [27][28][29] across the MPT, with a shift from 40 kyr obliquity-driven glacial cycles to 100 kyr paced glacial cycles [30][31][32] . The change in glacial-cycle length and magnitude across the MPT resulted in a change in sea-level variations, with sea-level fluctuation amplitudes increasing from 60-80 m before the MPT to ~120 m from MIS 22 (Fig. 3). The impact of the resultant scale of sea-level change on coastal systems across the MPT has typically been recorded only in raised marine terraces 23,33 . Our results from the Cooloola and Fraser Island coastal dune systems suggest a dramatic sediment-supply change driven by the switch to higher-magnitude sea-level variability at this time.

Sea level and sediment supply linked to dune building
The activation and stabilization of coastal dunes have been shown to be linked to sea-level variability and sediment supply [34][35][36][37] . Previous work in eastern Australia suggests a link between rising sea level and coastal dune formation due to the reworking of sediment stored on the continental shelf 20,22,38 . The presence of Late Pleistocene (~17.3 ka) parabolic dune forms at -60 m water depth, 40 km off the Fraser/Cooloola ( r e f .   39 also proposed that longshore transport was enhanced at the -60 m shoreline position due to a lack of bedrock headlands restricting northward transport in the region. Sediment stored in these mid-shelf settings was reworked shorewards during transgressions, resulting in large coastal dune-field formation during highstands 20,22,35,41 . North of ~33° S, the eastern Australia continental shelf becomes broader and shallower with lower slope, which allows depositional processes to dominate 35 . Due to the combination of shelf morphology, persistent northward longshore sediment transport and the relative tectonic stability of eastern Australia, a large seaward-thickening and prograding sediment wedge exists on the middle to outer shelf 35,[42][43][44] . This sediment wedge is estimated to be at least Miocene in age offshore of central New South Wales 42 . This mid-shelf sediment wedge would have been substantially disturbed during the MPT. Additional sea-level lowering by 40-60 m during periods of increased global ice volume would have exposed sediment that was previously stored in middle-to outer-shelf positions below wave base. This sediment was reworked shorewards 35,39 during periods of sea-level rise within 100 kyr glacial cycles. Sand that may have accumulated gradually over millions of years was transported into the inner continental shelf and near-shore coastal system during the first few 100 kyr cycles during the MPT. We interpret that the effect of this elevated sediment supply was to generate broad new coastal dune fields, specifically the Cooloola Sand Mass and Fraser Island, and probably the other sand islands of southeastern Queensland (North Stradbroke and Moreton Islands).
In this Article, we propose that shoreward shelf sediment reworking was enhanced as a direct result of the increased magnitude of sea-level change across the MPT. Translation of these previously stored mid-shelf sands landwards across the shelf led to increased sediment supply along the highstand coastline, which became available for deflation and aeolian transport as large transgressive dunes. Once this shelf sediment body was exhausted, the dune fields remained stable except for periods of parabolic reworking during the Late Pleistocene and Holocene highstands 22 .

The role of Fraser Island in shelf sediment transport and the GBR
The sand islands and dune fields that formed along the eastern Australia coastline during the MPT are not located randomly but occur where the continental shelf broadens as sediment travels northwards. Fraser Island is the northernmost and largest of the sand islands and acts as a 'groyne' that extends to the northeast across the continental shelf ( Figs. 1 and 4). It directs near-shore transported sand up to 80 km offshore and to within a few kilometres of the shelf break, where it is routed over the shelf and down submarine canyons to the abyssal plain 14 . This modern transit of sand from the end of Fraser Island to the continental slope is well documented 14 and is visible in satellite imagery (Fig. 4).
The southernmost emergent coral reef in southern Queensland is at Lady Elliot Island (24.11° S) (Figs. 1 and 4), about 80 km north of Fraser Island. However, this is not the southern ecological limit of corals on this coast; waters farther south are warm enough to support corals. Coral reefs are found in Moreton Bay (Fig. 1) at ~27.5° S, and there is evidence for somewhat more extensive mid-Holocene reefs 45 . These minor reefs occur in sheltered waterways behind barriers that block sand from entering, which would limit access to hard substrate and impact coral photosynthesis.
The absence of major reefs between Moreton Bay and Fraser Island is directly attributable to the high sand volume in the longshore drift system. We argue that before the formation of Fraser Island ~790 ± 90 ka (the oldest age obtained from this study), sand migrated north into the southern GBR and was retained on the inner shelf by the dominant southeast trade winds (Fig. 4). This is supported by Boyd et al. 14 , who showed that sediment travelling north of Fraser Island initially continues to the northwest before the strong ebb tidal currents from Hervey Bay entrain the sand and shift it back to the northeast towards Breaksea Spit (Fig. 4). Before Fraser Island developed, northward longshore transport would have interfered with coral reef development in the southern and central GBR. Distal sediment records from the south-central GBR (Ocean Drilling Program Site 1195; Fig. 1) support this interpretation, with a marked carbonate content increase and concomitant terrigenous and siliciclastic sediment decrease at ~700 ka (refs. 5,6 ). This change was interpreted to represent southern GBR initiation, with reef formation leading to inshore siliciclastic sediment trapping 5 . We propose, however, that it was the formation of Fraser Island that prevented longshore northward sand transport, which resulted in the observed decrease in terrestrial sediment transport to the south-central GBR thus allowing carbonate sedimentation to dominate. Fraser Island formation provided a necessary step for initiation of the southern and central provinces of the GBR during the MPT (Fig. 4). It is likely that this deflection of longshore sediment transport also contributed to reef development farther north along the GBR, but additional data are needed to test this linkage.
Several hypotheses have been proposed to explain GBR formation. These include (1) sea-surface warming due to the northward movement of the Indo-Australian plate 46-48 , (2) the prolonged and unusually warm MIS 11 interglacial (424-374 ka), which promoted global reef development 8,49,50 and (3) widening and deepening of continental shelves caused by longer-and higher-amplitude highstands across the MPT 8,51 . While northward movement of the Australian continent and subsequent transition of the Queensland coastline into tropical and subtropical waters was necessary for coral growth, this happened much earlier than the proposed development of the GBR 47,48 . Widespread reef development at MIS 11 is also not supported by more recent studies, which suggest an older GBR initiation by 450-670 ka (refs. 2,4,5 ). There is now mounting evidence for GBR initiation coinciding with both the MPT and formation of the southeast Queensland sand islands and dune fields.
Our results indicate that the shift in sea-level variability across the MPT led directly to the elevated coastal sediment supply that created Fraser Island and the Cooloola Sand Mass. Once emplaced, the orientation of Fraser Island influenced northwards sediment transport and provided a necessary precondition for coral development in the southern and central GBR. Fraser Island and the Cooloola Sand Mass represent a globally important archive of coastline geomorphic development in response to changes in Earth's orbital forcing during the MPT. This major coastal sedimentary system reorganization is not expected to be unique to southeastern Queensland, however, and should be examined in other passive-margin coastlines such as the Wilderness Embayment coastline of South Africa and the extensive coastline of eastern Brazil.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-022-01062-6.
Article https://doi.org/10.1038/s41561-022-01062-6 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Article https://doi.org/10.1038/s41561-022-01062-6

Drone image acquisition and orthophoto production
To determine the stratigraphy of the Rainbow Beach cliffs, high-resolution drone images were acquired to produce an orthophoto mosaic of the entire section. Images were acquired using a DJI Mavic Pro drone equipped with a 12-megapixel camera. A total of 211 photographs were taken covering the ~9 km length of the Rainbow Beach cliffs. Care was taken to ensure that the elevation and distance from the cliff remained relatively constant to reduce issues associated with stitching or distortion. This work was done only to produce a two-dimensional orthophoto of the cliff section; ground control points were not used.
Images were analysed and processed using the Agisoft Photoscan (http://www.agisoft.com) software to generate the point cloud and orthophoto mosaic of overlapping images following a standard workflow 57 and processing parameters. Orthophoto mosaics were exported as TIFF files and were viewed and analysed within ArcGIS Pro, and image annotation and editing were performed in Adobe Illustrator.

Field descriptions and sampling
To determine the stratigraphy and sedimentology of the Rainbow Beach cliffs, three sections (western, central and eastern) were logged in detail. At each site, sedimentological and pedological information was recorded, including visual grain size, sorting, roundness, bedding structures, mineralogy, Munsell colour and texture. Additional sub-samples were collected for further grain-size analysis. Pedostratigraphy was used to infer different dune activity phases. On the basis of podzolization, each pedological unit theoretically consists of E, B and C horizons (the C horizon is often absent). Each pedological unit can be grouped into an allostratigraphic unit, which represents a distinct period of dune deposition followed by stabilization and soil formation 58 . Each allostratigraphic unit is bounded by a contiguous discontinuity that represents a separate dune deposition period. At each location, the aim was to log the entire vertical section, although the steepness of the cliff face and presence of loose unstable sand made it impossible to reach all parts of the cliff. At all sites except for the Carlo Blowout section (Extended Data Figs. 1-4), only the bottom half or so of the cliff was logged at the face. However, with the aid of the orthophoto mosaics, the units identified at the Carlo Blowout can be traced and correlated laterally along the cliff 59 .

OSL sampling and analysis
Optically stimulated luminescence samples were collected from the exposed face of the major stratigraphic units at each site following standard procedures outlined by Nelson et al. 60 . Additional samples were collected away from the main sites to test the lateral continuity of allostratigraphic units. Each sample was collected by pounding a 20 cm steel pipe into the freshly cleaned outcrop surface. Representative dose rate samples were collected from the surrounding sediment within a 30 cm radius, and water-content samples were obtained in situ from the hole made by the sample tube.
Processing and analysis of the OSL samples was conducted at the Utah State University Luminescence Laboratory. Single-aliquot regenerative-dose 25,61 analysis was performed on small aliquots (~10 grains) of 180-250 µm quartz sand. Due to the intense leaching that has occurred in the E horizons, we chose small-aliquot dating instead of single-grain dating to reduce the influence of micro-dosimetry [62][63][64] . Samples were processed under dim amber (~590 nm) light conditions. Quartz sand (180-250 µm) was isolated by sieving and purified using 10% hydrochloric acid (HCl) to remove carbonates, peroxide to remove organics, heavy-liquid separation (2.7 g cm -3 sodium polytungstate) to remove heavy minerals and three 30 min treatments in 47% hydrofluoric (HF) acid to remove feldspars and etch the quartz. These steps were followed by a 37% HCl wash to remove any precipitated fluorites. Purity of the quartz grains was monitored using response to infrared stimulation and aliquots with signals greater than twice the background were rejected.
Dose rate calculations were performed on representative sub-samples using inductively coupled plasma (ICP) mass spectrometry (MS) and ICP atomic emission spectroscopy (AES) to determine the concentrations of K, Rb, Th and U in the bulk sediment surrounding the samples (Table 2), and internal dose rate 65 was calculated from the analysis of purified quartz sand (Extended Data Table 1). Moisture content for all samples was determined from the in situ sample collected from within the sample tube; however, due to drying effects on the sediment moisture content along the face of the exposed outcrops, an assumed value of 7 ± 2% was used to represent the average moisture history for all samples. Sediment radio-chemistry, cosmic-ray contribution and water-content attenuation were used to calculate dose rates following Guérin et al. 52 and Brennan et al. 53 . A cosmic-ray contribution was determined using the sample burial depth (from the top of the cliff) and the elevation and latitude and longitude of the sample site 54 .
Optical measurements were performed on small-aliquot (1 mm diameter, ~10 grains per disk) samples using Risø TL/OSL Model DA-20 readers with blue-green light-emitting diodes (LEDs) (470 ± 30 nm) as the stimulation source. The luminescence signal was measured through 7.5 mm ultraviolet filters (U-340) over 40-60 s (250 channels) at 125 °C with LEDs at 70-90% power (~45 mW cm -2 ) and was calculated by subtracting the average of the last 5 s (background signal) from the first 0.7 s (4 channels) of the signal-decay curve. Results of a preheat-plateau dose-recovery test 66 suggested that a 200 °C preheat for 10 s produces the best results for samples in this study. The luminescence signals decay rapidly and are dominated by the fast component of the signal 61 . For samples with <1 Gy equivalent dose (D e ), dose-response curves were fitted linearly between the zero dose and repeated regenerative doses. Where samples extended beyond the linear section of the regenerative curve, aliquots were fitted with either a saturating exponential or saturating exponential plus linear fit based on the fit to the dose-response curve for each sample to calculate individual aliquot equivalent dose values 25 . Example signal-decay curves and dose-response curves for the samples are presented in Extended Data Fig. 6.
Equivalent dose values were calculated using the Central Age Model of Galbraith and Roberts 55 using at least 16 accepted aliquots of quartz sand (see Extended Data Fig. 7 for radial plots of D e values). Aliquots were rejected if they had evidence of feldspar contamination, a recycling ratio beyond 20% of unity (<0.8 or >1.2), recuperation >1 Gy or natural D e greater than the highest regenerative dose given. Errors on D e are reported at two sigma standard error, and age estimates are reported at one-sigma standard error. Reported uncertainties include errors related to instrument calibration, dose rate and D e calculations, and errors were calculated in quadrature using the methods of Aitken and Alldred 67 and Guérin et al. 52 .

Palaeomagnetic sampling and analysis
The quality of palaeomagnetic recording is controlled by the magnetic domain state of magnetic particles, which depends on particle size, where uniformly magnetized single-domain particles are ideal recorders of magnetic signals, and larger multi-domain particles have poor palaeomagnetic recording capability. The particle size of sand exceeds the threshold size for multi-domain behaviour (several µm), which means that sands and sandstones are generally not ideal for palaeomagnetic investigations. By contrast, ferricretes contain abundant iron (oxy-hydr-)oxide nanoparticles that are more likely to be in the single-domain size range (tens to hundreds of nanometres depending on the mineral) and capable of carrying reliable palaeomagnetic signals. An important issue for palaeomagnetic studies of such weathering products is the timing of the iron mobilization events that created the ferricrete. If there is a relatively short time difference between the iron-leaching event and deposition of the host sand, palaeomagnetic Article https://doi.org/10.1038/s41561-022-01062-6 identification of geomagnetic polarity reversals can provide a useful minimum age for the host sand. In this setting, iron mobilization events are expected to happen quickly after dune deposition due to the warm subtropical climate, high annual rainfall and high permeability of the dune sands. Previous dating work by Ellerton et al. 68 and Walker et al. 20 shows that well-developed soil B horizons with prominent humic and sesquioxide accumulations (Bhs horizon) have formed in dunes that are less than 10 kyr in age.
Even when dune sands are coated by iron (oxy-hydr-)oxide pigments during iron remobilization, which could make them suitable for palaeomagnetic analysis, sand units at Rainbow Beach are friable, and samples usually did not survive transportation to the laboratory for palaeomagnetic analysis. Ferricretes tend to be sufficiently well cemented to provide coherent samples that survived transportation. Although ferricrete of this nature represents a different iron mobilization process from podzolization, these events are also not expected to post-date dune deposition by an appreciable amount. On the basis of the rates of soil formation in this region, we estimate that iron mobilization associated with ferricrete and the acquisition of remnant magnetization occurs 10-20 kyr after dune deposition. Palaeomagnetic samples were, therefore, taken from ferricretes in B or occasionally E horizons within sand units. Ferricretes were oriented by marking a vertical arrow on their outer surface and were semi-excavated to ensure that they were large enough to cut into multiple samples. Northward-directed arrows were marked on flat excavated surfaces of such ferricretes. Ferricrete samples were then removed from the outcrop and were placed in plastic bags for transportation to the laboratory. A few weakly cemented ferricretes did not survive transportation, but most remained intact. Samples were cut with a diamond-tipped saw blade in the laboratory and were placed into 2 cm × 2 cm × 2 cm plastic cubes for palaeomagnetic analysis. Markings were used to orient samples with respect to north; a 10.5° correction was made to compass readings to adjust for the present-day geomagnetic field declination at Rainbow Beach. Sample orientation was accurate only to ~5°. The purpose of this study is to identify the primary palaeomagnetic polarity of samples. Potential orientation errors of ~5° are small compared with the required precision.
Results of a pilot study on samples taken during initial palaeomagnetic sampling revealed that a dominant part of the palaeomagnetic signal of ferricrete samples from Rainbow Beach is carried by goethite (α-FeOOH), which has a maximum unblocking temperature of ~120 °C (ref. 69 ). The red colour of many ferricretes also indicates the presence of haematite (α-Fe 2 O 3 ), which has a maximum unblocking temperature of ~680 °C (ref. 70 ). Both minerals have high coercivities and are not expected to demagnetize substantially during progressive alternating-field (AF) demagnetization, which makes thermal demagnetization necessary to remove secondary overprints to identify the primary magnetization direction. Nevertheless, progressive AF demagnetization is not destructive and was carried out to assess whether contributions exist from lower-coercivity minerals. AF demagnetization was carried out at 23 steps of 2. remaining natural remanent magnetization (NRM) measured after each step using an inline demagnetization system attached to a 2-G Enterprises superconducting rock magnetometer system at the Black Mountain Palaeomagnetic Laboratory, Australian National University (ANU). As expected, the NRM was demagnetized only partially even in the highest applied AFs, so progressive stepwise thermal demagnetization was undertaken after completion of the preceding AF demagnetization treatment. The potentially dominant contribution from goethite means that multiple demagnetization steps should be applied below 120 °C. Most ovens used for thermal demagnetization are insufficiently calibrated to avoid thermal overshoot at such low temperatures. Thus, samples were analysed at the Geological Survey of Japan, where a Natsuhara-Giken oven with precise temperature control was used for thermal demagnetization with NRM measured on a 2-G Enterprises superconducting rock magnetometer system. Further samples were analysed at ANU after an ASC Scientific oven was modified to enable precise temperature control, with remanence measurements using a 2-G Enterprises superconducting rock magnetometer. Thermal demagnetization was carried out at steps of 80, 85, 90, 95, 100, 105, 110, 115 and 120 °C in Japan. Most samples were largely demagnetized at 120 °C, which confirms the importance of goethite in these samples. In other cases, the NRM remained partially demagnetized, so further demagnetization was carried out as needed at higher-temperature steps of 140, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 680 and 700 °C at ANU. For stably magnetized samples, a characteristic remanent magnetization (ChRM) direction was calculated using principal component analysis following the method of Heslop and Roberts 71 . For sites with multiple stably magnetized samples, site mean palaeomagnetic directions were calculated with errors associated with ChRM directions propagated into the mean directions following the method of Heslop and Roberts 72 .
Stable ChRM directions were obtained for ferricrete samples from Sites RB 41,42,43,45,47,50,51 and 52 (Extended Data Fig. 8, Extended Data Table 2 and Supplementary Figs. 1 and 2). Clear normal polarity results were obtained from sites RB 41-43 (lower dune package 2; Supplementary Fig. 1a-f) and sites RB 50-52 (dune package 3; Supplementary Fig. 2c-h). Clear reversed-polarity results were obtained from site RB 47 (dune package 2, 30 cm above the contact with dune package 1; Extended Data Fig. 8c-f). Noisy and less-clear results, including possibly transitional or reversed-polarity directions, were obtained from sites RB 44 and 45 (dune package 1; Extended Data Fig. 8a,b and Supplementary Fig. 1g,h). Results from site RB 49 are ambiguous, with some samples recording normal polarity (for example, Extended Data Fig. 8g,h), which could be a present-day field overprint, and other samples providing hints of an underlying reversed-polarity component (for example, Supplementary Fig. 2a,b). No final interpretation is provided here for site RB 49. The only samples with stable ChRM directions obtained from our pilot study were from ferricretes at Site RB 3, which yielded transitional ChRMs with reversed-polarity declinations (base of dune package 2; Extended Data Table 2).
Recording of exclusively normal polarity directions in the younger dune packages 2 and 3 is consistent with the OSL dates obtained for these packages and indicates that they were magnetized during the most recent Brunhes normal polarity chron (<773 ka (ref. 12 )). Recording of reversed polarity and transitional directions at several sites in lowermost dune package 2 and in dune package 1 is also consistent with OSL dates and suggests that the Matuyama/Brunhes palaeomagnetic reversal (773 ka) 12 , which provides an important time marker for the MPT, is recorded in this lower part of the dune sequence at Rainbow Beach.

Data availability
The datasets generated and/or analysed during the current study are available in the Bolin Centre for Climate Research Database repository at https://doi.org/10.17043/ellerton-2022-osl-queensland-1.
Article https://doi.org/10.1038/s41561-022-01062-6 Extended Data Fig. 8 | Example demagnetization results from this study that show a reversed, transitional and normal polarity for samples RB45, 47 & 49 from Rainbow Beach. Left-hand column: plots of NRM intensity with respect to demagnetization step with data points after alternating field demagnetization (AFD) in blue squares and thermal demagnetization (ThD) in yellow circles; green lines indicate data portions used for PCA fitting of ChRM directions. Righthand column: vector demagnetization diagrams from which declinations (blue circles) and inclinations (yellow squares) are determined by fitting a ChRM using PCA through multiple demagnetization points. Linear PCA fits are shown with orange lines for the inclination and green lines for the declination. Results shown above highlight example results for reversed, transitional and normal polarities for samples (a, b) RB45-5 (transitional polarity with southward declinations and downward inclinations); (c, d) RB47-10 (reversed polarity); (e, f) RB47-11 (reversed polarity); and (g, h) RB49-11 (a normal polarity sample from a site with generally unclear results). For remaining demagnetization results the reader is referred to the Supplementary Figs. 1-2.