Four boreholes ( Table 1 ) were drilled in the study area in September 2018; their coordinates and elevations were obtained using TOPCON RTK, Japan. Single tube rotary drilling was used, with a coring tube diameter of 100 mm. The core sections were sealed and transported back to the laboratory, where they were split longitudinally, photographed and the lithology described. One half was then sampled at 5-cm intervals and the other half was sealed and stored in a freezer. The locations of samples for OSL dating were determined based on the lithology. A total of 23 OSL samples and 996 samples for grain-size analysis were taken from the four cores.

OSL samples were taken by hammering in metal tubes vertically into the core face, immediately after splitting. Layers with a high silt content were selected for sampling. The tubes were packed in black plastic bags and immediately sealed with adhesive tape to prevent light exposure and moisture loss. Sample preparation was conducted under subdued red light. The unexposed material in the central part of the tubes was wet sieved to extract the 38–63 μm fraction, treated with 10% HCl and 30% H₂O₂ to remove carbonates and organic matter, respectively, and then etched with 35% H₂SiF₆ for two weeks to dissolve feldspars and then with 10% HCl to remove fluoride precipitates. The purity of the quartz grains was checked by infrared (IR, λ = 830 nm) stimulation. The isolated quartz grains were mounted as medium-sized (6 mm) aliquots in the centre of stainless steel discs using silicone oil.

The single-aliquot regenerative-dose (SAR) protocol (Murray and Wintle, 2000) was used for equivalent dose (D_e) determination. The luminescence was stimulated by blue LEDs (λ = 470 ± 20 nm) at 130°C for 40 s using a Risø TL/OSL-DA-20 reader with 90% diode power and was detected using a 7.5 mm–thick U-340 filter (detection window 275–390 nm) in front of the photomultiplier tube. Irradiations were carried out using a ⁹⁰Sr/⁹⁰Y beta source in the reader. The preheat plateau test, dose recovery test, recuperation ratio and recycling ratio were analyzed on representative samples to select a suitable preheat temperature and to check the suitability of the SAR protocol.

In the current study, the concentrations of U, Th and K were measured by neutron activation analysis. The water content (mass of moisture/dry mass) was determined by weighing the sample before and after drying, with ± 5% uncertainties. The cosmic-ray dose rate was estimated for each sample as a function of depth, altitude and geomagnetic latitude (Prescott and Hutton, 1994). The dose rates and OSL ages were calculated by Dose Rate and Age Calculator (Durcan et al., 2015) and are listed in Table 2 . In order to facilitate the comparison between sedimentary and historical records, the luminescence age (LA) in this paper can be converted into calendar age (CA) by the following method: CA = 2018 (AD) -1000 * LA (ka); If the calculation result of CA is greater than zero, the calendar age use the notations AD, otherwise use BC.

Grain-size distributions were measured using a Malvern Mastersizer 3000 laser particle-size analyzer in the institute of coastal research of Ludong University, China. The pre-treatment procedure consisted of the removal of organic matter and carbonates by the addition of 10% H₂O₂ and 10% HCl, respectively, followed by dispersal using 10 ml of 0.05 mol/L (NaPO₃)₆ and treatment in an ultrasonic vibrator for 10 min (Lu and An, 1997).

The mean grain size ranges from silt (4-63 μm) to very fine sand (63-125 μm), with silt the dominant component in most samples. The down core clay content (<4 μm) is relatively uniform, generally < 20%, while the sand content (63–125 μm) shows the largest range of variation (0–90%). The sand layers rich in shell debris were found in boreholes cores ZKH (8-9m), ZKA (5.6-6.5m) and ZKB (3.7-4.2m), and the grain-size results also show that these layers mainly comprise very fine sand ( Figures 2C–E ). The grain-size depth profiles of the four borehole cores generally show two stages of variation. The lower part is generally coarser with a higher sandy component than in the upper part, reflecting a significant change in the sedimentary environment. The coarser sediments in the lower parts of the cores are mainly gray and dark gray, with cross bedding and wavy bedding ( Figures 2H, I ); while the sediments in the upper parts are mainly yellow silt without bedding structures ( Figures 2A, B ) and compacted silty clay ( Figures 2F, G ). According to the grain size, color and sedimentary structures of the sediments, the lower gray sediments may be delta-front facies, which are usually grayish or blackish due to reduction environment under the water; while the yellowish silts in the upper parts are possibly delta-plain facies, which are usually yellowish or brownish due to the oxidation environment formed by frequent exposure to air.

The OSL dating results are listed in Table 2 . Most of the samples show similar luminescence characteristics and normal distributions for the De values, suggesting that the medium-grained quartz (38–63 μm) was well bleached prior to deposition, which is consistent with a detailed study of the applicability of the OSL dating method to coastal sediments of the southern Bohai sea (Chen et al., 2013). Previous studies have shown that ¹⁴C dating of coastal sediments is prone to age overestimation due to the carbon pool effect (Stanley and Chen, 2000); however, reliable late Holocene ages can be obtained using quartz OSL dating (Li et al., 2018). As shown in Figure 3 , the OSL ages of the four boreholes generally change according to the stratigraphic depth sequence, with the age range of 0.12–4.60 ka. Based on the OSL ages and depth changes in sedimentary facies, the stratigraphy of the four boreholes can be sub-divided and correlated, as shown in Figure 4 . The sequences can be divided into four sedimentary units according to the isochronic planes of the OSL dating and the major stages of delta development recorded by the historical documents (Gao et al., 1989): Holocene transgressive strata (Before 40 BC), the ancient delta (40 BC–893 AD), the abandoned delta (893 AD–1855 AD), and the modern delta (1855 AD to present). The sand layers within the cores, combined with the OSL chronology, could be used to determine the location of the coastline at various times as one of the evidences: in 40 BC (ZKH), 893 AD (ZKA), and 1855 AD (ZKB).

The spatiotemporal distribution of the former coastline of the YRD since 40 BC, when the Yellow River captured the Jishui River channel, has been reconstructed using various forms of evidence: sand layers containing shell fragments in three of the boreholes ( Figure 4 ), the distribution of saltworks during the historical period ( Figure 1A ), ancient shell beach ridges and buried shell beach ridges ( Figure 1A ), and the results of previous studies (Guo, 1980; Zhuang et al., 1991; Ji, 2009). The coastlines of 40 BC and 925 AD were determined based on topography, borehole data (including the sand layers with shell fragments in core ZKH), and the ¹⁴C ages of ancient shell ridges (Zhuang et al., 1991). The coastline of 893 AD was determined based on sand layers with shell fragments in core ZKA. The coastlines of 1855 AD and 1934 AD and the modern YRD were determined based on the results of Guo (1980), together with field investigations and the borehole core data (including the sand layers with shell fragments in core ZKB). The results show that the abandoned Yellow River delta was mainly distributed in the central plain of the Yellow River delta between the 1934 AD coastline and the saltworks pre-dating the Qin dynasty (Before 256 BC). As shown in Figure 1A , the four boreholes constitute a cross section from land to sea. The cross section is parallel to the left bank of the Daqing River, which is roughly located along the axis of the shallow-buried abandoned YRD and spans the three coastlines of 40 BC, 893 AD, and 1855 AD. The spatiotemporal distribution of sediment deposition recorded by the boreholes reveals that the shallow-buried abandoned delta has the following three geographic features:

(1) Inverse terrain gradient from land to sea

As shown in Figure 4 , in 893 AD, when the delta was abandoned, the elevation of the terrain of the ancient delta decreased from southwest to northeast; and in 1855 AD, before its burial, elevation of the delta decreased and then increased from southwest to northeast. The abandoned delta deposits abruptly become thicker from southwest to northeast, showing a slender horn-shape. The OSL age of the highlands on the seaward side of the delta is 1019–1799 AD (1–0.22 ka), and that of the main body is 1379–1799 AD (0.64–0.22 ka). Therefore, different from the pattern of the terrain of the modern YRD, which has an elevated long axis (containing the river channel) and lower sides, with higher elevations upstream and lower elevations downstream. It was also different from the ancient delta of 893 AD that dipped gradually from land to sea, with the elevation of the terrain of the abandoned delta high on both sides and low in the middle, along the northeast-southwest direction, with height differences of ~3.3 m and ~1.0 m, respectively. The slopes of the southwest and northeast parts are 2.03 ‱ and 1.44 ‱, which are respectively ~3 times and ~2 times greater than the slope of the modern deltaic plain (0.75 ‱) ( Figure 4 ). The northeast side is a low and gently sloping coastal highland distributed along the 1855 AD coastline, which is located between the 925 AD and 893 AD coastlines, with a width of ~20 km. The top of the highland is ~1.01 m higher than the lowest part of the depression, and ~2.59 m higher than the top of the 893 AD beach sediment layer in the ZKA borehole, forming an inverse gradient from the plain depression to the coastal highland in the SW–NE direction ( Figures 1 , 4 ).

Previous studies have shown that in the five counties of the YRD during the Ming and Qing dynasties, Lijin and Guangrao, which were along the long axis of the delta, experienced seven and five severe storm surges, respectively, during 1374 AD–1855 AD. This was much lower than at Wudi (16) on the west side, and at Zhanhua (11) and Shouguang (17) on the east side (Wang, 2010). However, the landward part of the coastal highland was a sea salt–producing area for more than 500 years (1289–1855 AD) in the middle and late stages of the abandoned delta ( Figure 1 ). Before the 16th century, salt in the YRD was produced using topsoil with a high salt content and shallow underground brine as raw materials, rather than directly using seawater (Ji, 2009). This indicates that before 1855 AD there was a large-scale positive terrain barrier along the coast of the abandoned delta and a shallow-buried high concentration of underground brine in the landward plain.

(2) Sandy coast and deep estuarine channel

The 40 BC and 1855 AD coastlines in the central plain of the YRD demarcate the land-sea boundary before the ancient delta formed and before the abandoned delta was buried, respectively ( Figure 1 ). However, there are two shorelines corresponding to the ancient delta when it was abandoned. The age of the first is 893 AD, while the ¹⁴C age of the lower shell ridges of the second is 925 ± 70 AD (Chen, 1982), with the intervening distance of ~20 km ( Figure 1 ). Considering the possible error in the ¹⁴C dating results and that the ¹⁴C age of the other site of the younger shoreline is 835 ± 90 AD (Chen, 1982), it can be determined that the geomorphic ages of the two coastlines are the same. Therefore, at the beginning of the formation of the abandoned delta there was a compound coast, and offshore shell ridges were developed 20 km away from the continental coastline ( Figure 5C ). A low-energy sedimentary environment of coastal shallow lakes and marshes was developed between the ridges and the land, which subsequently evolved into the above mentioned low elevation and gently sloping highland, and a single sandy coast developed outside the highland.

After the Yellow River captured the Jishui River and entered the sea in 40 BC, the middle and lower reaches of the ancient Jishui River gradually evolved into the Daqing River, entering the sea along the southern depression of the ancient delta. After 1077 AD, the lower reach of the Daqing shifted northeastward towards the sea via Lijin, incising the abandoned delta and forming a new estuary ( Figure 1 ). The estuary of the Daqing River is 15-20 m deep with the abandoned delta plain as its first terrace, which is 6-9 m higher than the low-water level. The channel is a single rectangular section 30 m wide and 6-9 m deep, with the tidal boundary being above Lijin (Gao et al., 1989). Given the earliest establishment of Tiemenguan wharf in the Jin Dynasty (1115-1234 AD), shortly after 1128 AD, it can be inferred that the abandoned delta was incised after the diversion of the Daqing River. Subsequently, a single, meandering, rectangular-shaped deep estuary gradually formed. Until 1855 AD, it was an important port for salt transport, but shortly after 1855 AD it became silted up, abandoned and eventually buried.

(3) Depth of sediment accumulation

In addition to the incised landform of the estuary channel of the Daqinghe River, the abandoned delta plain is an accumulation geomorphic body stacked on top of the ancient delta. However, the sediment types, thickness and the sequence of the different geomorphic components vary substantially ( Figure 4 ). The thickness of the coastal highland after abandonment exceeds 4.0 m. The protruding part is 1.5-2.0 m–thick with abundant shell fragments and thin silt layers belonging to lake–marsh facies. The silt is black or dark gray in color with a high organic matter content; it changes gradually upward to grayish yellow continental silty clay or clayey silty sand ( Figure 3 ). The OSL chronology for the upper part of the ZKB and ZKA boreholes shows that the sedimentary ages at the base and top of the highland are 1019 ± 100 AD and 1799 ± 30 AD, respectively. The former age is later than the formation of the abandoned delta, and the latter age is earlier than the burial of the abandoned delta. The protruding part of the highland recorded in the ZKB borehole was formed during 1379 ± 60 to 1799 ± 30 AD, in the late stage of delta abandonment.

The abandoned delta in the plain depression accumulated a ~0.5-1.0 m–thickness of continental silty clay sediments with a sediment deposition rate less than that of the northeastern coastal highland, and also less than that of the overlying modern delta and the underlying ancient delta. According to the chronology of the ZKH borehole core, on the west side, the average deposition rate of the plain depression is only 0.87 mm/yr; the rate for the modern delta is 25.37 mm/yr, that for the ancient delta is 5.31 mm/yr, and that for the highland ZKB borehole core during the same period is 4.38 mm/yr ( Figure 3 ). The sediments for almost the whole of the ZKH borehole core and the lower part of the ZKE borehole core in this layer are compact, with brownish-red and brownish-yellow colors and a low sand content ( Figures 2F, G ). Additionally, they contain numerous black carbon fragments or horizontal carbon bands, and a small amount of Fe-Mn rust spots, termed the “laterite layer” by local villagers during well digging. This indicates that the vertical aggradation of the abandoned delta plain was slow and pedogenesis relatively strong.

The period between the emergence of the abandoned delta and its burial was ~1000 years (893AD–1855 AD), spanning the interval from the Medieval Warm Period to the Little Ice Age, when there were major changes in geomorphic processes, including a reduction in the sediment supply, and changes in storm surge intensity. These changes are discussed below.

(1) Reduction in sediment supply

Delta construction is the result of the input, transportation and accumulation of river sediment, and changes in the river estuary will inevitably have a profound effect on the geomorphic evolution of an abandoned delta (Peng et al., 2010). After 893 AD the Yellow River changed its course from the ancient YRD, entering the Bohai Sea to the north of the abandoned delta (Guo, 1980). After 1128 AD most of the runoff and sediment flowed southward, capturing the channel of the Huai River to the South Yellow Sea. After 1495 AD all of the water and sediment load of the Yellow River entered the South Yellow Sea via the Huai River; previously, only a small amount of water and sediment entered the sea via Lijin. During this period, at ~1077 AD, the Daqing River (It originated from the Shandong province) flowed into the Bohai Sea through the abandoned delta via Lijin, forming a new estuary (Chen, 1982). Therefore, during the intervals of 893–1077 AD, 1077–1128 AD, 1128–1495 AD and 1495–1855 AD, the runoff and sediment supply to the abandoned delta underwent major changes. The influence of the Yellow River on the abandoned delta gradually weakened, and after 1495 AD it ceased completely, leaving the abandoned delta affected solely by the erosive estuary of the Daqing River, which has a much lower sediment load (Chen, 1982). Therefore, the existence of shell ridges and sand layers with shell fragments in the borehole cores indicates that, with the reduction of sediment supply from the Yellow River, the effects of wave erosion on the abandoned YRD increased, transforming the ancient delta from a muddy coast to a sandy coast.

(2) Changes in climate and storm surges

Changes in climate, storm surges and wave activity substantially affect coastal environments and hence the geomorphic evolution of coastal estuaries and deltas (Yang et al., 2019). Speleothem δ¹⁸O and δ¹³C records from Kaiyuan cave, 150 km southwest of Lijin, provide a regional climatic record during 892–1894 AD (Wang et al., 2015). The 590 years before 1482 AD comprised the Medieval Warm Period (MWP) when the summer monsoons was strong and stable, and the 412 years after 1482 AD comprised the Little Ice Age (LIA) when the summer monsoon was weak and unstable ( Figure 6 ). Strong storm surges prevailed around the Bohai Sea during the past 2000 years (Feng, 1982; Liu, 1989). Based on the observing hydrodynamic and sediment data, the simulation results of an abandoned Yellow River Delta lobe, revealed that local storms enhanced tidal residual currents and considerably weakened tidal shear front, intensified the sediment resuspension and dispersal (Fan et al., 2020). Under storm conditions, surface suspended sediment concentration was 2 times as large as that in normal conditions. Moreover, the weakened tidal shear front was in favour of the sediment dispersal in the river mouth and southward diffuses to the central area of the Laizhou Bay (Fu et al., 2021). Laizhou Bay, where the modern YRD is located, experiences a storm surge once every three years on average during the recent 50 years; during intervals of strong wind and wave action, the storm surge height increases by 2.11–5.74 m and lasts for 24-72 hr (Wang et al., 2016), providing energy and an accumulation space for coastal landform development. According to statistics from historical documents (Liu and Zhang, 1991), there were nine storm surge disasters in the ancient delta area of the Bohai Sea during the interval of 933 years (40BC-893AD), but as many as 34 storm surges affected the abandoned delta during the interval of 962 years (893AD-1855AD). During the 962 years of the abandoned delta, for the 589 years prior to 1482 AD there were only 3 storm surges, while in the 373 years after 1482 AD there were as many as 31 ( Figure 6 ). As shown in Figure 1 , the 1855 AD coastline is located on the landward side of the 925 ± 70 AD coastline, indicating that the abandoned delta had migrated inland in the late abandonment stage, forming a coastal highland. Therefore, the frequent and intense storm surge activity after 1482 AD, during the Little Ice Age, may have promoted the formation of the coastal highland–depression system of the abandoned delta.

Based on the foregoing analysis of the processes affecting the evolution of the YRD, combined with the sedimentary evidence from the borehole cores and the results of previous studies (Guo, 1980; Zhuang et al., 1991), the geomorphic evolution of the shallow-buried abandoned YRD during 893–1855 AD ( Figure 5 ) was reconstructed. These evolutionary stages are summarized below.

In brief, after the reduction in sediment supply in 893 AD, the simplex muddy coast was transformed to a compound sandy coast. The coast did not retreat significantly, but the topographic relief increased and the slope became steep. The large-scale ancient deltas, which first developed some 2000 years ago, were rapidly abandoned some 1000 years ago. They were then exposed during the interval from the Medieval Warm Period to the Little Ice Age and then rapidly buried to a shallow depth. This process was accompanied by significant changes in estuary types, which were unique amongst river deltas on a global scale. Their geomorphic characteristics were completely different from those of the modern Yellow River and the abandoned YRD in northern Jiangsu, and they also differed from those of the Luan River delta in the northern Bohai Sea (Xu et al., 2020) and the ancient estuary lagoon coast in the southern part of Laizhou Bay (Han et al., 2002). They differed from the classical sandbar-lagoon system but also differed from the abandoned (sub-) deltas of other major rivers.