Abstract
The world’s mega-deltas are extremely important from a human perspective and attract considerable effort to reveal their evolution, growth-related driving forces, and human impacts. Here, we report a case study on the Holocene deltaic evolution of the Yellow River, through development of a conceptual model, which is compared with paleo-proxy to analyze the forcing acting on the delta. The main conclusion is that superlobe switching was modulated by the 1500-year cycle. Cooling in Mongolia in response to strong Bond IRD events, which is coincident with warming in eastern China due to a strong Kuroshio Current, enhances the meridional temperature gradient, which then increases cyclone frequency and activates dust storms and terrestrial erosion throughout the catchment. Enhanced erosion supplies great amounts of material to the Yellow River and causes channel evulsion and superlobe development, expressed as dominant 1500-year cycle. At the same time, summer monsoon and solar forcing are uncorrected with deltaic evolution on these timescales. Therefore, we conclude that Holocene dynamics of the delta on a millennial timescale was dominated by winter cyclone activity across northern China and Mongolia.
Keywords
Introduction
Mega-deltas, such as those of the Nile, Mississippi, Mekong, Yangtze (Changjiang), and Yellow River (Huanghe), are important from a human perspective because they are home to major cities housing millions of people (e.g. Nicholls, 1995; Woodroffe et al., 2006), are important economic centers (e.g. Niou, 2002), and encompass diverse ecosystems (e.g. Macintosh, 2005; Sanlaville and Prieur, 2005). The complexity of these systems makes them vulnerable to environmental variability (e.g. Nicholls, 2004; Woodroffe et al., 2006). Sea level rise, human activity, and climate change are the three major factors affecting the vulnerability of these mega-deltas (Intergovernmental Panel on Climate Change (IPCC), 2007). It is thus important to understand the linkages between naturally occurring deltaic processes, anthropogenic impacts, and the internal and external driving forces that combine to generate this variability (e.g. Nicholls and Cazenave, 2010; Syvitski et al., 2009; Wang et al., 2011).
The Yellow River delta is characterized by extremely large sediment loads and contributes significantly to the globally integrated, riverine sediment supply (Milliman and Syvitski, 1992). The Yellow River delta has attracted considerable attention from researchers in various fields, such as sedimentary evolution (e.g. Cheng and Xue, 1991; Gao et al., 1989; Liu et al., 2009; Qiao et al., 2011), deltaic geomorphology (e.g. Saito et al., 2000; Xue, 1993), and modern sedimentary process (e.g. Chu et al., 2006; Li et al., 1998; Wang et al., 2011; Wright et al., 1986). Due to this extensive research base, our understanding of the evolutionary history of the Yellow River delta and associated superlobe formation has been gradually unfolding.
Superlobe formation occurs in response to a shift in the river’s terminal channel when high sediment loads in the lower reaches (Saito et al., 2000), choke off the old channel, leading to evulsion and the creation of a new terminal channel. The initiation of superlobe formation in the Yellow River dates to the mid-Holocene stabilization of sea level around 6000 14C yr BP (Hori and Saito, 2007; Stanley and Warne, 1994; Xue, 1993). A total of 10 paleo-superlobes of the Yellow River (Figure 1) have developed (Xue, 1993). They have been identified by shell ridge studies around the western coast of the Bohai and Yellow Sea (Xue, 1993; Xue and Cheng, 1989). Three of these superlobes (numbers 1, 6, and 7) have been confirmed by stratigraphic and sedimentary studies (Saito et al., 2000; Xue, 1993). Nine of the superlobes formed on the western coast of the Bohai Sea, while the remaining one was laid down on the western coast of the Yellow Sea between AD 1128 and 1855 (Xue, 1993), when the Yellow River fed into the Yangtze River (Liu et al., 2010). Human impacts on deltaic processes due to changing land use in the past 2000 years and later the construction of major river dams have only been significant recently (Milliman et al., 1987; Ren and Zhu, 1994; Saito et al., 2001). The present-day Yellow River delta began forming in AD 1855 after a large-scale engineering project shifted the course of the Yellow River from the Yellow Sea back to the Bohai Sea. This current Yellow River delta is referred to as superlobe 10 (Figure 1).

(a) Climate systems probably affecting the northern China. The Yellow River delta and its catchment are labeled as a triangle and oval, respectively. IRD: ice-rafting debris, AMO: Atlantic meridional overturning, KC: Kuroshio Current, and YSWC: Yellow Sea warm current. Southeast and southwest monsoons are two systems in summer. (b, c, d) Study area, paleo-superlobes (Saito et al., 2000; Xue, 1993; Xue and Cheng, 1989), paleo-shorelines, and YDZ-3 location. Ten paleo-superlobes of the Yellow River delta in the Holocene were labeled as numbers. Modern marine surface sediments in the Yellow River estuary (d1; Ren et al., 2012) and Xiaoqinghe estuary (d2; Yi et al., 2012) were collected for analysis.
Of the three major factors influencing deltaic development and vulnerability (sea level rise, human activity, and climate change), the impact of climate change on deltas is the least well understood. This lack of understanding arises from a combination of poor sample spatial coverage, inadequate geochronology, and poorly quantified proxies for deltaic evolution (e.g. Marriner et al., 2012).
In this study, we investigate sedimentary records extracted from a new borehole located in the southern Yellow River delta, focusing on how climatic change affected the natural process of deltaic evolution during the Holocene. Radiocarbon dating is combined with deltaic proxy studies for this purpose. This work enhances our understanding of the natural background processes operating on mega-deltas and their potential impacts on climatic change and earth landform evolution.
Methods
Study area
The Bohai Sea is a semi-enclosed marginal sea on the northeast coast of China (Figure 1), with an average water depth of 18 m. It is connected to the northern Yellow Sea by the narrow Bohai Strait, which exceeds 60 m depth (Marine Geology Laboratory, Institute of Oceanology, Chinese Academy of Sciences (IOCAS), 1985). The dominant hydrodynamic factors in the Bohai Sea are tidal currents, waves, and ocean currents (Wei et al., 2004). The Bohai Sea was likely exposed during the Last Glacial Maximum (18,000–26,000 cal. yr BP), when sea level was 120–130 m lower than present (Hori et al., 2001; Yang and Lin, 1993). Regional sea level rose up to its present level by ~8000 cal. yr BP and to its most recent highstand 6000–7000 cal. yr BP, after which sea level oscillated slightly before falling to the present level (Xu, 1994; Zhao et al., 1979).
The principal source of sediment to the Bohai Sea is the Yellow River. The Yellow River, which originates in the northeastern part of the Tibetan Plateau, flows across the Chinese Loess Plateau, before discharging into the western Bohai Sea. The sediment load carried by the Yellow River has increased dramatically during the last 2000 years as a consequence of cultivation and de-forestation on the Chinese Loess Plateau (Milliman et al., 1987; Ren and Zhu, 1994; Saito et al., 2001), and the annual quantity of sediment is about 1.08 × 109 tons (Milliman and Syvitski, 1992). However, in more recent times, human activity has had a profound impact on the sediment supply from the Yellow River to the Bohai Sea (Wang et al., 2007, 2011). Moreover, it is estimated that with a fluvial discharge rate of 7.67 × 109 m3 yr−1 and a sediment load of 0.278 × 109 ton yr−1 (Xu, 2002), the Yellow River delta stands at the threshold between erosion and aggradation.
During the early Holocene transgression around the Bohai Sea, the sedimentary environments adjacent to the drilling site (see below) changed from rivers and lakes to tidal flats and estuaries (Marine Geology Laboratory, IOCAS, 1985; Xue et al., 1995), after which the paleo-environment shifted to a deltaic one around ~6000 14C yr BP (Marine Geology Laboratory, IOCAS, 1985; Xue et al., 1995; Yu et al., 1999; Zhao et al., 1978).
YDZ-3 core
The YDZ-3 core, located on the southern margin of the modern Yellow River delta (37°27′N, 118°55′E; Figure 1), was drilled in 2009 by East China Normal University. Sediment was transported along the northern margins of the paleo-superlobes (Figure 1c) to the drill site, producing an expanded section. Prior to AD 1855, the paleo-shorelines was landward of the drill site (Figure 1d), although the superlobe prograded past the drill site by AD 1934. This single location is well situated to monitor the Holocene development of the Yellow River delta. Because of its unique location, any single site further to the northwest would not suffice for this purpose since it would be situated upstream of several of the paleo-superlobes.
The length of core obtained from the drill site was 34.7 m, and the recovery rate was 91%. An exceptional section of the record with 98% recovery raised from 18.0 to 7.7 m below the core top was selected here for study. This 10.3-m-thick interval is undisturbed and well constrained by radiocarbon dating. Because no sedimentological evidence was observed for erosional unconformities in this interval, we divided it into two units based on lithology: (1) the deeper unit, spanning 18.0–15.4 m below the core top, is composed of a basal layer of dark gray and black organic-rich clay (20 cm thickness) which is overlaid by grayish silt and clayey silt, with mollusk debris. This lower unit is interpreted as a shallow-sea system. (2) The shallower unit from 15.4 to 7.7 m below the core top is composed of yellowish fine silt, which is similar in composition to the modern sediment in the Yellow River delta. This unit is interpreted as a deltaic environment. The vertical variation of these sedimentary facies is depicted in Figure 2.

Stratigraphic columns, radiocarbon dates, grain size spectrum, and color changes in CIELAB space of the borehole YDZ-3. L* describes the lightness between black (0) and white (100), while a* and b* denote the red–green and yellow–blue chromaticity, respectively.
Radiocarbon dating and age model
Seven peat samples from the core were extracted for radiocarbon dating (Figure 2). All accelerator mass spectrometry (AMS) radiocarbon dating measurements were conducted at the AMS Laboratory of Peking University in the China in 2010. To remove reservoir effects and stable carbon isotope influence, radiocarbon dates were calibrated with the Calib 6.1 program (Stuiver et al., 2009) using ΔR value of 143 ± 59 years (a regional reservoir data for the Yellow River delta and Laizhou Bay; Wang and Fan, 2005; Wang et al., 2004; Table S1, available online).
The core site was located on the southern margin of the modern Yellow River delta and the western coastal plain of Laizhou Bay. The Yellow River, as well as local rivers such as the Zimaigou and the Xiaoqinghe, transports sediment into the Bohai Sea (Figure 1d). When the sediment from the Yellow River decreased due to superlobe shift, local rivers would become the main source at the core site. Because of this unique geomorphological setting and the lack of any significant depositional hiatuses observed during the core sampling, we applied a second-order polynomial fit, which considered the 95.4% confidence level for each radiocarbon date during the fitting process (Yi et al., 2013), to develop an age model for the depth interval from 18.0 to 7.7 m in the YDZ-3 core.
Sediment grain size
A total of 207 grain size samples were measured at a 5-cm sampling interval. The grain size samples were pretreated with 10–20 mL of 30% H2O2 to remove organic matter, washed with 10% HCl to remove carbonates, rinsed with deionized water, and then placed in an ultrasonic bath for several minutes to facilitate dispersion. The grain size spectra of the remaining terrigenous material were measured using a Malvern Mastersizer 2000 laser-particle size analyzer at East China Normal University. A total of 100 grain size classes between 0.3 and 300 µm were exported for further analysis.
Sediment grain size is a powerful proxy for paleoenvironmental reconstruction because depositional interpretation varies with sedimentary grain size and composition. Because grain size spectra represent mixtures of sediment delivered by multiple processes, to identify the processes controlling grain size variation, varimax-rotated, principal component analysis (VPCA) is often employed. This method partitions the variance in the grain size dataset into sediment input components that can be interpreted in terms of processes (e.g. Darby et al., 2009; Yi et al., 2012).
We seek to determine whether surface sediment and down core samples are derived from the sample population. As an initial step to test this underlying assumption, we performed separate VPCAs on marine surface and suspended sediment datasets (Ren et al., 2012), on marine surface sediment from the Xiaoqinghe estuary (Yi et al., 2012), and on the core data presented in this study. All of these datasets were later combined for a single, unified analysis. Loess grain size collected from the Xifeng Profile on the Chinese Loess Plateau (Hao et al., 2008) was also compared for reference. This allows us to test the null hypothesis that different processes control sedimentation in various datasets (Darby et al., 2009; Yi et al., 2012). If this were the case, the components extracted from the different datasets would not match.
Mineralogy by diffuse spectral reflectance
Here, we employ the CIELAB color space in which any sediment’s color can be expressed by L*a*b* values. While the CIELAB indices provide useful stratigraphic information (e.g. Ortiz et al., 2009), they integrate information over too broad a spectral range to be diagnostically linked to specific mineralogy in a way that is theoretical as opposed to empirical.
Measurements of diffuse spectral reflectance (DSR) were collected at 5-cm resolution in the Sedimentology Lab at Kent State University, USA, from the wet, sampled, and well-preserved sediment using a Minolta CM-2600d ultraviolet/visible (UV/VIS) spectrophotometer (wavelength range: 400–700 nm; resolution: 10 nm; spot size: 3 mm) following the methods of Ortiz et al. (2004). The sediment sample was covered/wrapped in a single layer of Gladwrap™ plastic wrap to prevent contamination of the instrument’s integration sphere during the measurement process. We used Gladwrap for this purpose, as is standard practice in the Ocean Drilling Program and Integrated Ocean Drilling Program, because it has been shown to absorb less light on the UV end of the spectrum than other commercially available wraps (Balsam and Deaton, 1991).
We employ VPCA to identify the contribution of different mineral assemblages to the down core color variations. The resulting dataset was decomposed by VPCA and analyzed by visible derivative spectroscopy following the methods of Ortiz (2011). The VPCA was calculated using the correlation matrix of the center-weighted derivatives of the sediment DSR (Ortiz, 2011; Ortiz et al., 1999, 2004) from the samples collected from the YDZ-3 core.
Mineralogy by x-ray diffraction
As an independent check on our visible derivative spectroscopy-based mineral identifications, we compare the down core components extracted from the DSR measurements on the YDZ-3 core with standardized mineral assemblages estimated by x-ray diffraction (XRD) from representative samples in the same core. A total of 12 down core samples from the YDZ-3 core were prepared for XRD analysis according to the methods described by Liu et al. (2004).
The XRD clay mineral study was carried out on the <2-µm fraction, which was separated by conventional Stokes’ settling after the removal of carbonate and organic matter by acetic acid (15%) and hydrogen peroxide (10%), respectively. Clay minerals were then identified by XRD using an X’Pert PRO, PANalytical XRD instrument (40 kV and 40 mA) at the Institute of Geology and Geophysics, Chinese Academy of Sciences. Two XRD runs were performed for each sample, following air-drying and ethylene–glycol salvation at 60°C for 12 h. Identification of clay minerals was made according to the position of the (001) series of basal reflections observed on the XRD diagrams (Moore and Reynolds, 1997).
Results and analyses
Chronology of the YDZ-3 core
The radiocarbon dates from the YDZ-3 core increase stratigraphically with depth, indicating that within the chronological resolution of the depth–age pairs, the sediment accumulation rate follows a consistent relationship throughout the Holocene. The second-order polynomial fit, displayed in Figure 2, and the individual regression parameters are all significant at p < 0.01, demonstrating this age model is statistically significant and robust. Other estimates of Holocene sediment thickness around the Yellow River delta are up to 20 m (e.g. Liu et al., 2009), consistent with the value of ~18 m that we estimated for the YDZ-3 core.
Grain size analysis
Variations in the clay, sand, and silt fractions are subtle in the older sedimentary unit but more pronounced in the younger stratigraphic unit (Figure 2). The grain size VPCA results of the individual datasets (Figure 3) indicate that the leading components, GSC-1 and GSC-2, account for roughly 75% of the variance in the suspended sediment, marine surface sediment, and down core sediment samples in each dataset. While there are differences in the minor components (GSC-3, GSC-4, and GSC-5) between the datasets, the principal difference is often the relative ranking of these components, rather than the addition of a distinct component. Because there is evidence that the same five components are present in each dataset, we conclude that the various datasets have essentially the same data structure and can be analyzed together for comparison and interpretation.

VPCA results from different types of samples: (a) YDZ-3 sediment (this study), (b, c) estuary surface and suspended sediments (Ren et al., 2012), (d) Xiaoqinghe marine surface sediment (Yi et al., 2012), and (e) all of the sediments. (f) Loess samples from Xifeng Profile (Hao et al., 2008) were also displayed for reference. The percentages labeled in plots represent variances of each component from VPCA.
While the dominant components of grain size have not changed, their relative percentages have varied through time. To evaluate the environmental information included in these stable components for further analyses, we calculate the weighted average of the two leading components to form a new series, GSC-12, which captures 74.9% of the overall grain size variability:
The weighting coefficients are simply the percentage of variance associated with each component. This approach filters out the minor constituents from the grain size spectra using the data adaptive filters defined by the first and second components (Ortiz, 2011; Yi et al., 2012). High values of GSC-12 indicate greater abundance of fine grains supplied by increasing Yellow River sediment input, while lower values of GSC-12 represent a coarsening of the marine sediment when Yellow River sediment influx decreases.
In general, GSC-12 variation can be grouped into two stages: low variability before ~6200 cal. yr BP and high variability after (see Figure 8a). From 9400 to 6200 cal. yr BP, GSC-12 is characterized by low variation without any distinct changes. During the early Holocene, GSC-12 has an average value of −39.4 with a standard deviation of 14.6. In contrast, from 6200 to 1600 cal. yr BP, GSC-12 exhibits greater variation, with an average value of −15.1 and standard deviation of 56.5. This zone of high variability can be partitioned into four sub-stages: 6200–4300, 4300–3400, 3400–2500, and 2500–1600 cal. yr BP. Moreover, all of the transitions between each sub-stage are abrupt. For example, GSC-12 is only −86.1 at 4249 cal. yr BP but becomes 40.8 at 4208 cal. yr BP, after which it sharply decreases from 78.4 at 4009 cal. yr BP and then to −100.6 at 3893 cal. yr BP.
Mineralogy by DSR and XRD
Proxy of mineralogy
Sediment brightness (L*) oscillates between values of 24.88 and 50.97 between 15.4 and 18.0 m depth and has a prominent valley (indicating black) at 17.7–17.9 m correlating to a 20-cm-thick, organic-rich layer. The brightness then increases into the upper portion of the sediment record with much less variation. The blue–yellow contrast (b*) and the red–green contrast (a*) exhibit features similar to L*, but without the valley at 17.7–17.9 m.
The results show that the three leading rotated DSR components account for 90.0% of the variance in the correlation matrix of the reflectance derivative spectra. The leading component accounts for 33.4% of the variance, while the second and third components account from 31.7% and 24.9% of the variance, respectively.
To identify the minerals or mineral assemblages associated with each component, we used Pearson’s correlation coefficient to calculate the linear correlation between each component and the center-weighted, first derivative spectra of known minerals and mineral mixtures measured in the laboratory at Kent State University or available online from version 5 of the USGS Digital Spectral Library (Figure 4). The first principal component (DSC-1) correlates to a mixture of chlorite + smectite (r = −0.97), the second (DSC-2) to a mixture of illite + goethite (r = 0.92), and the third (DSC-3) to glauconite + actinolite (r = 0.95). The first derivative spectrum for pure goethite has principal peaks at 535 and 435 nm (Deaton and Balsam, 1991), which are similarly observed in Figure 4b. Similarly, the other minerals identified by visible derivative spectroscopy have peaks consistent with those found in the standards to which they are compared.

Varimax-rotated factor loadings and DSR first derivates of mineral standards plotted as a function of wavelength. (a) The first principal component (DSC-1) correlates to a mixture of chlorite + smectite (r = −0.97), (b) the second (DSC-2) to a mixture of illite + goethite (r = 0.92), and (c) the third (DSC-3) to glauconite + actinolite (r = 0.95).
As an independent check on the visible derivative spectroscopy analysis, additional results were obtained by XRD analysis from 12 representative sediment samples (Figure 5). The identified minerals are smectite (4.6–4.9 Å), chlorite (6.1–6.3 and 18.7–18.9 Å), illite (8.8–8.9 and 17.6–17.8 Å), mixed layers of kaolinite–chlorite (12.2–12.5 and 24.8–25.3 Å), quartz (20.8–29.0 Å), and mixed layer of illite–quartz (26.5–27.0 Å). The result of the XRD-based clay mineral analysis is consistent with the DSR results in this study and with XRD and DSR results for loess deposits (e.g. Bronger and Heinkele, 1990; Liu, 1985).

Variations of (a, b, c) DSR components with (d) 12 representative samples for XRD analysis. Position of each XRD sample was marked as solid circles and empty ovals. The identified minerals are labeled.
Down core variation of mineralogy
Because the down core time series are derived by VPCA, the three leading series are not intercorrelated, but rather represent an orthogonal decomposition of sedimentary variations (Figure 5). This simplifies interpretation of the results because each component includes only information that is independent from the others in the VPCA result, but which can be related to specific mineral assemblages reflecting distinct provenance or processes.
Chlorite + smectite content, expressed as DSC-1, shows pattern of change similar to GSC-12 throughout the Holocene. From 9400 to 6200 cal. yr BP, the chlorite + smectite content is characterized by subtle variation. In contrast, from 6200 to 1600 cal. yr BP, chlorite + smectite content varies considerably with four sub-stages identified at 6200–4300, 4300–3400, 3400–2500, and 2500–1600 cal. yr BP. As with the GSC-12 record, all of the alternations between sub-stages in the DSC-1 record are abrupt and coincident in their timing.
Illite + goethite content, expressed as the DSC-2 series, is initially relatively low, but increased sharply after 6200 cal. yr BP. Cyclic changes then dominate the illite + goethite record, which exhibit three peaks around 6140, 4630, and 3130 cal. yr BP. Transitions between high and low contents are all abrupt, occurring within ~50 years.
Glauconite + actinolite content, expressed as the DSC-3 series, shows a two-stage variation: high fluctuation before ~6200 cal. yr BP and little variability after, which is in opposition to the GSC-12, DSC-1, and DSC-2 records. From 9400 to 6200 cal. yr BP, glauconite + actinolite content increased rapidly from a relatively low value, with two peaks at 8480 and 7590 cal. yr BP, respectively. From 7590 cal. yr BP, the content decreased stepwise until 6200 cal. yr BP. In younger portions of the record, the glauconite + actinolite content varied with much higher frequency than observed in the other records.
Spectral analysis
Spectral analysis in the frequency domain was conducted to identify major cyclicities present in all four of the independent series (GSC-12, DSC-1, DSC-2, and DSC-3) for the portion of the Holocene that was recovered. We employ the ARAND software (Howell et al., 2006) to conduct the spectral analysis of these proxies. The results indicate that the 1500-year cycle is the dominant periodicity in all the series except for DSC-3 (Figure 6). In addition to the statistically significant millennial and longer cycles, there were shorter centennial-scale cycles present, in the range from 100 to 400 years.

Spectral analysis of four series, that is, (a) GSC-12, (b) DSC-1, (c) DSC-2, and (d) DSC-3, by the multi-taper method, which was from the ARAND software (Howell et al., 2006). The 1500-year and its half cycles (600–800 years) were highlighted as gray bars. Significant minor components were also labeled. The horizontal axis is in natural log scale.
Discussion
Mineral composition
Chlorite, smectite, and illite are common clays in the sediment of the Bohai Sea and Chinese Loess Plateau (e.g. Bronger and Heinkele, 1990; Liu, 1985). Goethite often forms through the weathering of other iron-rich minerals (Klein et al., 1993). Furthermore, goethite is widely found in sediment from the Chinese Loess Plateau (e.g. Heller and Evans, 1995; Liu et al., 2006) and the Chinese marginal seas (e.g. Chen, 1989). Considering the composition of loess deposits and large amount of sediment input by the Yellow River, it is plausible that these clays were mainly derived from the Chinese Loess Plateau. Glauconite and actinolite are both greenish minerals, but indicative of different sedimentary processes. Glauconite is considered a diagnostic mineral, indicative of a marine depositional environment on continental shelves (Klein et al., 1993) and is often identified in sediment from the Chinese marginal seas (e.g. Chen, 1989). Actinolite, however, is derived mainly from metamorphic rocks (Klein et al., 1993) and is reported in high concentration in Yellow River sediment (e.g. Lin et al., 2003; Yang et al., 2009), indicating a genetic relationship to the Yellow River discharge. Therefore, the mineralogical interpretations of the leading three components derived by derivative spectroscopy can be interpreted and tied to sediment material that is regionally common, indicating that the DSC-1 and DSC-2 are genetically correlated to the Yellow River and the DSC-3 might reflect the competition between marine influence and terrestrial inputs in the Holocene, related to the authigenic formation of glauconite in the marine environment and the hydrodynamic redistribution of dense, fluvial-derived actinolite by wave and current action on the shelf.
Correlating proxy to paleo-superlobes of the Yellow River
Regional sea level rose up to the present level at ~8000 cal. yr BP and to the highest level at 7000–6000 cal. yr BP, after which, it oscillated slightly before falling to the present level (Xu, 1994; Zhao et al., 1979). Regional tectonism during the Holocene was calm without significant uplift or subsidence in the study area (Marine Geology Laboratory, IOCAS, 1985). From these observations, we hypothesize that the hydrodynamics associated with deltaic erosion was much less than sediment input during the Holocene, resulting in ample accommodation space. This hypothesis is consistent with the measured progradation of the delta prior to construction of major river dams (Milliman et al., 1987; Ren and Zhu, 1994; Saito et al., 2001). Furthermore, the YDZ-3 core was located downstream relative to each paleo-superlobe with the exception of PSL-9 and PSL-10 (Figure 1c), and the paleo-shorelines prior to AD 1855 were considerably landward (Figure 1d), providing ample room for progradation of the southern margin of the paleo-superlobes toward the study site. Because the tidal and geostrophic flow pathways in the Bohai Sea are generally cyclonic, sediment from the paleo-superlobes of the Yellow River would have been carried toward the study site. Due to its distal position, the YDZ-3 core site is thus well situated to record the major changes in the deltaic evolution of the Yellow River during the Holocene. Any single location further to the northwest in the Bohai Sea or along its coast would not suffice for this purpose because it would be upstream of several of the paleo-superlobes.
Empirically, estuaries and deltas develop at river mouths in response to transgressive and regressive phases, respectively (Boyd et al., 1992), modulated by variations in sediment accumulation rate (e.g. Hori and Saito, 2007). Deltas typically are classified according to the interplay of their sediment supply and their primary control on deposition: river-, wave-, or tide-dominated systems (Galloway, 1975). Most natural systems are intermediate between these three end-member states. Although deltaic dynamics are complicated, all actively prograding deltas have at least one common attribute (Coleman, 1982; Wright, 1977): a river supplies clastic sediment to the receiving basin more rapidly than marine processes can remove it. Accordingly, we describe this process as follows:
Here, Δ represents deltaic changes; Ss and Se indicate sediment supplied to and eroded from the depositional basin, respectively; t and a represent the temporal and spatial variabilities of the deltaic processes, respectively; and ε represents stochastic factors, including those not explicitly defined in the conceptual model. The model implies that the probability of delta formation varies as a function of sediment supply and erosion, both of which vary with time and space, with some random component added.
For a given delta in which river input is much larger than erosion, such as the Nile, Mississippi, or Yellow River deltas, the probability of delta existence should be positively related to the riverine sediment supply:
Here, X indicates the matrix of environmental proxy values related to the riverine input, P(Δ) represents the probability function for deltaic formation, and C is a constant coefficient representing sediment erosion due to hydrodynamics. Now, this conceptual model states that the greater the river input, the larger the probability of deltaic development, modulated by sediment erosion of course. We can thus reconstruct the probability of deltaic evolution of the Yellow River delta on various spatial and temporal scales depending upon river inputs.
As discussed in the previous section, three of the four proxies can be employed to indicate river input from the Yellow River: (1) GSC-12, derived from the sediment grain size spectra, indicates variations in the abundance of fine-grained sediment supplied by the Yellow River and deposited in the delta; (2) DSC-1 and DSC-2, derived from sediment mineralogy, represent the riverine input of the Yellow River derived from different parts of the watershed; and (3) the fourth proxy, DSC-3, might represent the competition between marine influence and terrestrial inputs and their redistribution in response to wave and current action.
According to Eq. 3, the greater the river input, the larger the probability of deltaic development, modulated by sediment erosion. We thus take GSC-12, DSC-1, and DSC-2 as proxies for net sediment input in Eq. 3 and the probability of deltaic formation, P(Δ). We assume that C is minimal during most of the Holocene because of the ample accommodation space and predominant sediment flow pathways. Comparing the sub-stages identified from the four time series indicative of Yellow River input with the known timing of paleo-superlobe formation on the Yellow River delta, we obtain a reasonably consistent fit between the two datasets over the time span from 6200 to 1600 cal. yr BP (Figure 7).

Correlating proxy to paleo-superlobe variation of the Yellow River. (a) GSC-12, (b) DSC-1, (c) DSC-2, (d) DSC-3, and (e) 10 paleo-superlobes of the Yellow River delta. Seven of them, labeled as PSL-1 to -7, are confirmed in this study and highlighted in color bars.
Can this conceptual model explain the paleo-superlobe variations of the Yellow River during the Holocene? Due to the post-glacial transgression, sea water first entered the Bohai Sea around 11,600 cal. yr BP (Liu et al., 2007), and its western coastline receded rapidly, reaching a position 50–90 km west of the present coastline during the maximum transgression at 6000–7000 cal. yr BP (Wang et al., 1986; Zhao, 1986; Zhao et al., 1978). From ~9400 to 6200 cal. yr BP, GSC-12, DSC-1, and DSC-2 all exhibit low values, and in contrast, DSC-3 reaches essentially high values by ~8500 yr BP, indicating marine dominance in sediment hydrodynamics at this area at that time. The Yellow River delta began forming on the western shore of the Bohai Sea at ~6000 14C yr BP (Saito et al., 2000; Xue, 1993), when sea level rise slowed (Hori and Saito, 2007; Stanley and Warne, 1994). In this study, GSC-12, DSC-1, and DSC-2 all sharply ascend but DSC-3 decreases at ~6200 cal. yr BP, providing independent evidence that the Holocene delta of the Yellow River probably initiated at that time. Following the deltaic initiation, GSC-12, DSC-1, and DSC-2 maintained high values and then gradually declined, until ~5000 cal. yr BP, when the superlobe shift from the Lijin (PSL-1) to the Huanghua (PSL-2) occurred. Similar agreements between superlobe shifts and proxy variations correlate with the published timing for PSL-3, PSL-5, PSL-6, and PSL-7 as reported by Xue (1993), confirming the timing of these events during the Holocene.
Processes during the formation of the Huanghua (PSL-2) and the Shajinzi (PSL-4) superlobes differed from the others as inferred from their unique sediment composition and texture. During these periods, while illite + goethite content increased, fine grains and chlorite + smectite content remained very low without any significant variations. Our results provide additional documentation for the timing of these events, which were not well constrained by prior studies (Saito et al., 2000; Xue, 1993) but raise additional questions as to why sedimentary processes were different at these times. Additional evidence will be required to further constrain the factors that contributed to the differences observed for PSL-2 and PSL-4 relative to other paleo-superlobes.
In summary, the sedimentary proxies GSC-12, DSC-1, and DSC-2 can be employed to indicate Yellow River deltaic evolution during the Holocene. While the first two DSC components relate to the Yellow River provenance, the third component, DSC-3, details the competition between marine influence and terrestrial inputs, which separates it from the discussion that connects GSC-12 and DSC-1 and DSC-2 to the conceptual model of delta growth. According to these three proxies, the deltaic evolution in the Holocene displayed a dominant 1500-year periodicity with an associated 600–800 half-cycle (Figure 6), inferring a common forcing for the various components on larger spatial scales, which will be discussed later.
Driving forces of Holocene deltaic dynamics
The Yellow River is the second greatest river in the world in terms of sediment load and annually transported 1.08 × 109 tons of sediment to the sea prior to damming (Milliman and Syvitski, 1992). It is estimated that the Yellow River added ~2300 km3 of sediments to the sea during the last ~2000 years (Milliman et al., 1987). Such a vast amount of sediment eroded from the basin has made sediment supply the predominant process in the Yellow River delta formation. Moreover, the proxy analysis present here indicates that although the sediment is mainly from the Chinese Loess Plateau, after transportation and redeposition, the marine sediment grain size structure differs significantly from that of the Chinese Loess Plateau, integrating and recording new information arising from other depositional processes. Therefore, we postulate a direct relationship between the deltaic dynamics and the average erosion in the catchment, with stochastic processes in the system determining which part of the delta grows or erodes.
The catchment of the Yellow River is affected by several major climatic systems (Figure 1a). The East Asian summer monsoon brings moisture into the Yellow River catchment from the south, enhancing sediment delivery. The Kuroshio Current (KC) modulates the winter climate of eastern China, including the lower reaches of the Yellow River. In addition, through the Westerlies and winter monsoon, the signal of North Atlantic and Arctic climate anomalies (as recorded by ice-rafting debris events or Bond IRD events) could be transmitted into the Yellow River catchment. Each of these climate-related processes could subsequently influence deltaic evolution of the Yellow River, and the IRD and KC are also known to be dominated by 1500-year oscillations within the Holocene (e.g. Bond et al., 2001; Darby et al., 2012; Isono et al., 2009; Jian et al., 2000; Sorrel et al., 2012).
Impacts from Bond IRD events and KC
Bond IRD events are expressed in the North Atlantic when large numbers of icebergs are discharged from the Labrador Sea and Baffin Bay. These cooling events weaken North Atlantic deep water formation and Atlantic Meridional Overturning through increased surface freshening, which enhances stratification and increases sea ice cover (Bond et al., 1997, 2001). A far-field effect associated with the Bond IRD events is that they strengthen the Asian winter monsoon circulation by influencing the northern Westerlies, increasing dust storms in northern China (e.g. Sun et al., 2011).
At mid-latitudes, marine processes within the continental marginal seas of the Northwest Pacific can play a crucial role in environmental changes (Liu, 2009). The KC, which originates from the North Equatorial Current in the western Pacific, carries warm, saline water to the north (Barkley, 1970) and affects winter climate in eastern China (Hsueh et al., 1997; Li, 1993). During winters with strong KC flow, the current would warm eastern China, with the opposite occurring during years with weaker KC flow (e.g. Hsueh et al., 1997; Jian et al., 2000; Li, 1993).
These two dominant influences likely integrate in a consistent way: cooling in Mongolia by strong Bond IRD events (e.g. Sun et al., 2011) is associated with warming in eastern China by strong KC flow. This enhances the meridional temperature gradient, which then increases winter cyclone frequency and dust storms in northern China (Qian et al., 2002). This teleconnection is consistent with the observation that Bond IRD cooling events in the North Atlantic are associated with warm, low-latitude sea surface temperature records (Darby et al., 2012), likely inducing storm activity in the North Atlantic (Sorrel et al., 2012), and which indicates that the linkages are hemispheric in nature. In this situation, wind-driven erosion in the catchment will be induced as observed in Pleistocene and modern Chinese loess studies (e.g. Zhu et al., 2004), and more materials would thus be delivered for deltaic initiation and development. However, because of the rapid sea level rise in the early Holocene, the Yellow River delta could not develop until the middle Holocene, when sea level rise slowed (Hori and Saito, 2007; Stanley and Warne, 1994).
Furthermore, it is noteworthy that the transition from one superlobe to another occurs when both IRD and KCI values suddenly increased from local minima, following a short stochastically influenced, phase offset between the deltaic and climatic proxies (Figure 8). For example, IRD decreased from 16.2% at 5250 cal. yr BP, reached a local minimum of 3.9% at 4900 cal. yr BP, and increased again to 9.3% at 4550 cal. yr BP; meanwhile, KC increased from a local minimum at 5220 cal. yr BP to a local maximum at 4912 cal. yr BP. During this period from minimum to maximum, PSL-1 shifted to PSL-2 around 4900–5000 cal. yr BP. Similar transitions can be found at PSL-2 to -3, PSL-3 to -4, PSL-4 to -5, and likely PSL-5 to -6.

Deriving forces of Holocene dynamics of the Yellow River delta. (a, b, c) Three proxies indicating deltaic evolution of the Yellow River developed in this study. (d) IRD, stack ice-rafting debris data from Bond et al. (2001). (e) Stack KC, integrated indicator from four time series by Jian et al. (2000), based on VPCA procedures done in this study. Arrows in (d) and (e) indicate possible temporal consistency between (d and e) the climatic variations and (a–c) shifts of the Yellow River paleo-superlobes (see details in text). Seven paleo-superlobes of the Yellow River, labeled as PSL-1 to -7, are highlighted in color bars. (f) Total solar irradiance (TSI) during the Holocene (Steinhilber et al., 2009) with the low-frequency component (>100 years). (g) Monthly sediment loads of the Yellow River at Lijin gauging station (Figure 1c) from 1950 to 2011 (data from the Yellow River Conservancy Commission, http://www.yellowriver.gov.cn/).
Consequently, impacts of climatic changes on superlobe shifts in the Holocene may be depicted as follows: relative warming in Mongolia by weak Bond IRD events, which is associated with cooling in eastern China by weak KC flow, reduces the meridional temperature gradient, which reduces cyclone frequency and dust storms in northern China. As a result, erosion in the catchment is reduced, causing sediment loads of the Yellow River to decrease following delta lobe switching. When the sediment loads in the Yellow River descend low enough, specifically to 0.278 × 109 ton per year, the modern critical threshold for Yellow River delta formation as estimated by Xu (2002), conditions in the delta lie close to the silting threshold. Below this threshold value, the river channels become self-adapting, interacting with floods, storm surges, or variations in tidal regime or wave action through enhanced channel erosion and greater potential for avulsion, as discussed by Fagherazzi (2008). Once established, if sediment loads suddenly increase, the equilibrium is disturbed, causing the channel to become overloaded with a greater potential for channel avulsion. New distributaries form in this way enabling superlobe migration.
Impacts from the Asian monsoon
Winter climate in Mongolia and northern China, which significantly mobilizes dust storms in China (Qian et al., 2002), is not the only potential factor affecting deltaic evolution of the Yellow River. The massive amount of precipitation associated with the summer monsoon plays an important role in sediment delivery and erosion.
The Yellow River flows across the Chinese Loess Plateau, where eolian deposits are easily eroded. Monsoonal climate is one of the most influential factors controlling this erosion (Liu, 1985). When the winter monsoon strengthens, regional vegetation retreats and deserts activate (e.g. An et al., 1990, 2001; Guo et al., 2002; Liu, 1985), enhancing erosion within the Yellow River catchment. There are only eight time slices (850, 465, 100, 7, 0.7, 0.3, 0.15 kyr BP, and the present) reporting erosion rate across the Chinese Loess Plateau (Zhu et al., 2004). These data are not sufficient to assess whether erosion rate provide the dominant control on Yellow River delta evolution. Research also establishes that 1500-year periodicities were relatively weak in the Asian monsoon during glacial periods (e.g. Sun et al., 2010, 2011) and even much weaker during the Holocene (e.g. Wang et al., 2005). Although six alternations between loess and soil development in the Holocene were reported in the upper reach of the Yellow River (Chen et al., 1999), there is no temporal consistency between those alternations and the timing of paleo-superlobe shifts measured in this study. Additionally, we compared our proxies of the Yellow River delta with various monsoon indicators, but little consistency was observed at a millennial timescale (Figure S1, available online).
Chinese loess is deposited during the winter and experiences pedogenetic process during the summer (Liu, 1985). Because of this anti-phase relationship between summer and winter climates in the Asian monsoon on millennial temporal scales (e.g. Xiong et al., 1996), it is possible that those phase differences within the Asian monsoon cause the 1500-year component to be muted or hidden by alternations between loess and soil, potentially resulting in the weak relationship observed between the Asian monsoon and deltaic evolution.
Solar influence on deltaic evolution
There have been a number of prior studies supporting the hypothesis that variations in solar irradiance are linked to millennial-scale climate variations, such as the Bond IRD events (e.g. Bond et al., 2001), Asian monsoon (e.g. Wang et al., 2005), and KC (e.g. Jian et al., 2000). When comparing total solar irradiance (TSI; Steinhilber et al., 2009) with our proxies for the deltaic evolution of the Yellow River, no direct correlation was observed (Figure 8). The Yellow River lobe phases PSL-2, -3, and -4 appear to be correlated with relatively high irradiance of the sun, but other lobe phases did not show a similar relationship. Similarly, we did not find any persistent relationship between either the Bond IRD events or the KC record when they were compared with solar activity (Figure 8). This observation is consistent with new Holocene results related to North Atlantic storm dynamics (Sorrel et al., 2012) and the Arctic Oscillation (Darby et al., 2012), in which no persistent phase association between those climate processes and solar activity was observed on millennial timescale.
In the case of the Yellow River deltaic system, three processes might contribute to this complex relationship: (1) although erosion mainly occurs during winter, most of the eroded materials are carried by rivers during the monsoonal season due to the higher water discharge (Figure 8g). (2) Summer precipitation across the catchment could also enhance erosion during flood periods. In addition, considering that the ENSO system acts as a mediator for solar influence on the climate system’s low-latitude heat engine (Clement et al., 2001; Emile-Geay et al., 2007; Marchitto et al., 2010) and that there is a significant connection between ENSO and rainfall in northern China rainfall (e.g. Wang, 2006), which influences the Yellow River discharge, there should be some relationship between ENSO and deltaic dynamics of the Yellow River. Whether the teleconnection to ENSO is direct or indirect, and how it relates to solar forcing, remains a topic that needs to be explored.
Conclusion
We describe how the climatic system influenced landform process in the Yellow River delta during the Holocene. Our main conclusions are as follows:
Abrupt transitions in grain size, illite + goethite content and chlorite + smectite content in the YDZ-3 core all agree well with the timing of Holocene paleo-superlobe migration of the Yellow River, expressing a dominant 1500-year cycle.
Deltaic migration on the Yellow River is correlated with the paleo-climatological records of winter cyclone activity: cooling in Mongolia in response to strong Bond IRD events is coincident with warming in eastern China due to a strong KC, which enhances the meridional temperature gradient and increases cyclone frequency, which activates dust storms and terrestrial erosion throughout the catchment. Enhanced erosion supplies greater amounts of material to the Yellow River, causing channel evulsion and superlobe development. However, little consistency was observed from summer monsoon and solar forcing variations, likely suggesting other linkages may exist.
Therefore, we confirm that the Yellow River delta initiated at ~6200 cal. yr BP. Holocene dynamics of the delta on a millennial timescale was dominated by winter cyclone activity across northern China and Mongolia.
Footnotes
Acknowledgements
The authors wish to thank Drs Chunxia Zhang, Wenxia Han, LiangCheng Tan, Xin Zhou, and Hao Long for discussion. LY, SC, and JDO designed the study; SC collected samples; LY, GC, JP, FL, and YC conducted experiments; and LY, SC, JO, and CD wrote the paper. All authors contributed to data interpretation and provided significant input to the final manuscript.
Funding
This research was supported by the National Natural Science Foundation of China (41402153, 40925012), National Basic Research Project of China (grants: 2010CB951202, 2012CB821900), and the China Postdoctoral Science Foundation (grant: 2013T60164). Part of this work was completed while the lead author conducted a post-doc at Kent State University.
