Holocene Evolution of the Pearl River Delta: Mapping Integral Isobaths and Delta Progradation

1. Introduction

River deltas and estuaries are dynamic environments in which marine and fluvial processes come together to form complex and transient morphologies. Their internal sedimentary architecture, formed in environments in which fluvial, estuarine, and shallow marine depositional processes have fluctuated over both space and time in response to allogenic controls (e.g., sea-level rise, tectonics, climate change), autogenic processes (e.g., channel avulsion) [1,2,3,4,5] and human activities [6], makes up an important part of the geological record. Many of the world’s largest deltas are densely populated and are becoming increasingly vulnerable to the influences of natural (e.g., climate and sea level) and non-natural (e.g., human activity) factors. While infrastructural works such as coastal protection schemes and the development of underground real estate is ongoing in delta plains, it is particularly important to thoroughly understand the composition of Quaternary sedimentary formations.

The Pearl River Delta (PRD) is the largest delta along the coastline of southeastern China and has been the subject of a great deal of research into deltaic formation and evolution over more than a century. The study of the sedimentary evolution of the PRD plain began with the impressive work of Huang et al. in 1982 [7]. Since then, numerous studies of its sedimentology, geomorphology, paleoclimate, and pre-historic human activity have been carried out [8,9,10,11,12,13,14]. However, most of these have not been able to use core data at a high spatial resolution to reconstruct the paleo river courses and incised valleys that formed during the glacial period. In recent years, more and more detailed studies have been concerned with the formation of PRD in response to postglacial sea level rise [14] and human impact on the changing shoreline during the late Holocene [6,15,16,17]. In particular, proxy data from multiple cores in the PRD show a rapid transformation from marine facies to alluvial facies during the late Holocene which may, in part, indicate the impacts of human activity [18]. However, more comprehensive studies at a larger spatial scale are needed to support the inference that human activities, which have intensified since 3-2 ka, led to accelerated development of the PRD plain [6,15,19].

Many of East and Southeast Asia's largest population and commercial centers are located in coastal regions on rivers and deltas. Both natural processes, such as river discharge, sea-level rise, and human activities have had an ever-increasing impact on river deltas [20]. The river network in the PRD is one of the most complex in the world [21]. Existing research suggests that the past river system with radial and beaded patterns formed during the late Pleistocene, and the river network of the delta today basically follows such pattern [22,23]. To build on this understanding of deltaic plain geomorphology, river channels, and coastline shifts in the PRD, a synthetic study of the spatial distribution of Holocene deposits, based on a new dataset at a higher spatial resolution of overall deltaic areas, is needed.

For this study we established a dataset of over 2800 sediment cores from the PRD, and used a GIS-based 3D geological model incorporating newly-acquired digital borehole data to create isobath maps of all Quaternary and Holocene deposits in the PRD. A series of reconstructed paleo-coastline maps revealed the delta’s evolution and progradation since the maximum marine transgression. Comparison between the present-day river courses and the reconstructed paleo channels of incised valleys formed during the LGM allows for a more detailed study of natural and human-induced coastline shifts and land growth.

2. Regional setting

The PRD, located at the northern boundary of the South China Sea, is the second largest delta in China after the Yangtze River Delta. Developed within a basin confined by the dominant NW-SE-trending normal faults, the PRD intersects a set of older NE-SW-trending normal faults. Compared to other deltas along the Southeast Asian coast [24,25,26], Quaternary sediments in the PRD area are relatively thin, and their formation was influenced by two major marine-terrestrial sedimentary cycles that have occurred since the late Pleistocene. As a result, the two marine sediment units separated by a weathered clay or fluvial sand are commonly found in the study area [27]. The Quaternary sediments in the study region consist mainly of deltaic deposits composed by of marine-terrestrial sediment sequence, and alluvial deposits from the lower reaches of the rivers which were laid down before the delta plain formed (Figure 1). This paper deals primarily with the deltaic areas formed by Quaternary marine-continental sedimentary cycles (including the maximum transgression range during the Holocene), and excludes the areas of downstream pure alluvial and peat deposits from multiple tributaries near the delta margin.

The PRD is surrounded by hills and mountains, and encompasses numerous paleo-river mouth islands. River water input and sediment discharges come from three main streams: the West River (WR), North River (NR), and East River (ER), of which the WR is the largest with a drainage area of some 350,000 km². Other smaller rivers include the Tan River (TR) and Liuxi River (LXR). All of these tributaries converge in the delta area to form a very complex delta-plain network of rivers [28]. Average river network density is as high as 0.8km/km2 within the delta plain. The three main rivers (i.e. WR, NR, and ER) branch out across the delta plain into over 324 small distributaries and streams close to the land-ocean interface, amounting to a total channel length of about 1600 km [29]. Most watercourses in the delta flow into the estuaries of Lingding Bay and Modaomen Bay. The principal tributary river outlets at the front of the delta are Humen, Jiaomen, Hongqimen and Hengmen to the northeast, and Modaomen, Jitimen, Hutiaomen and Yamen to the southwest (Figure 1b).

3. Data and Methods

3.1. Borehole data

To illustrate the spatial distribution of the Quaternary deposits and basement topography in the PRD, we collected ca. 20,000 pieces of borehole data from multiple sources, including data newly-generated for this study, the existing literature, and unpublished documents from engineering drilling projects, to build a borehole database. All data were examined, and any cores whose lithological interpretation was insecure due to imprecise or ambiguous logging were excluded from the database. Reasons for excluding data included lack of coordinate information, insufficient lithological description, and insufficient drilling depth. Following this selection process, 2819 pieces of borehole data were included in the database and creation of a digital geological model (Figure 2). The database includes two types of datasets. The first of these (538 cores) contain detailed lithological descriptions, multiple proxy analyses data, radiocarbon dating results, and some high-resolution core photographs. In general, the data in this dataset were logged at a fixed interval resolution, for texture of sediments (clastics and organics), color, visible macrofossils such as brackish and fresh water mollusks or shell fragments (for common species, with an indication of the genus), and plant debris (e.g. charcoal, wood). Data from the unpublished new cores drilled with the support of land and resources protection and governance of Guangdong Province (No. 2017201) and the detailed information found in the literature from the past 40 years were included in this dataset. The second dataset includes a substantial volume of filtered borehole data, most of which were provided by Guangdong Geological Survey Institute and various engineering projects, and generally contains location coordinates, principal litho-stratigraphic descriptions and, in many cases, hydrological engineering test results.

3.2. Analysis of sediment cores and stratigraphic synthesis

For the 538 boreholes conducted for this study and collected from the literature, the data included not only all available lithological, stratigraphic and fossil/microfossil information, but also chrono-stratigraphic data (in particular dating results from 14C, OSL, and ESR). 13 new sediment cores were drilled from the delta plain to enable detailed analysis of their sediments, geochemistry, and microfossil content. This analysis was complemented by reference to representative and well-dated records in the published data. Synthesising the findings from this dataset enabled us to refine our analysis of the spatiotemporal distribution of Holocene sediments. A thorough revision of all the available stratigraphic logs and geotechnical tests improved the existing map of the thickness of Quaternary deltaic sediments. Overall, the dataset established during this study has been proven useful in determining the major sedimentary units of the late Quaternary.

To address the issue of complicated spatial facies changes, we examined all information from the 2819 cores in the dataset and summarized the stratigraphic sequences into six major sedimentary formations within Quaternary deposits as a whole, within which the Holocene comprises three major sedimentary facies formations (further details are given in the results section). In this study, only the results of mapping the Holocene sediment thickness are shown.

3.3. Method of isobath mapping

The present-day topographic map of the study area was produced using data derived from the USGS/NASA hole-filled seamless SRTM data at a 90×90 m resolution [30]. Areas of differing elevation were colored in ArcGIS (ESRI, version: 10.3) using the ‘symbology’ function (Grey color ramp, use hill-shade).

Sediment thickness maps are a fundamental tool in Quaternary geology. They represent spatial thickness variations and thickness trends of a given age. The borehole data discussed above formed the basis for the modelling procedure which generated the maps. All of the available data were stored and analyzed in a GIS environment using ordinary block-kriging Kriging interpolation technique through the software ArcMap [31,32]. The interpolation procedure was carried out to obtain the best surface fitting the borehole data, and to adapt for each individual stratigraphic unit. The results derived from variograms and search-neighborhoods (both key elements in the interpolation) were combined specifically. Elevation contour lines were drawn using the same software with the raster extraction algorithm ‘Contour’ (interval set to 5 m). By overlaying the mask layer (Pearl River Delta Plain region) with the bounding surface depths of different strata recorded in the boreholes, we were able to obtain the stratigraphic isobath map of each sedimentary unit. Given the scarcity of borehole data in the sea (particularly in Lingding Bay), the isobath map only shows the land portion of the delta plain (Figure 2).

In order to create a four-dimensional landscape model for the Quaternary evolution of the PRD, the changing spatial extents of the main sedimentary facies and the evolution history of the delta plain had to be established. For this paper, we first identified the depth of landscape surfaces of the weathering bedrock, allowing us to map the 'pre-delta' topography upon which the Quaternary delta was laid down (Figure 2b). Then, the late Pleistocene and Holocene bounding surfaces (i.e. spatial mapping of the Holocene deposits) were established based on their sedimentological identification in boreholes and the geological data contained within the chronostratigraphic database (Figure 2c). The spatial extents of the different environments within the delta were mapped through an in-depth analysis of the chronostratigraphic data which took into consideration the sedimentary units within which each dating result was hosted.

4. Results

4.1. Holocene stratigraphy and lithofacies changes

Tectonic, climatic, and sea-level changes are widely accepted as important controlling factors in the development of the Pearl River delta plains [33]. Quaternary deposits in the study area have previously been summarized and grouped into different stratigraphic units according to the regional heterogeneity of sedimentary facies [6,8,27,34,35,36,37]. The spatial variety of sedimentary facies and the complexity of the data sources mean that the stratigraphic divisions for the Quaternary sediment sequence varies between areas and literatures. In general, seven sediment formations have been identified (Table 1), which can be grouped into two major marine-terrigenous sedimentary cycles which occurred during the late Pleistocene and the Holocene, respectively. Based on sediment and lithofacies changes in the Holocene we are able to distinguish sequences of marine transgression and regression, two important processes in deltaic evolution which reveal how the coastline has moved over time. The Holocene lithostratigraphy and inferred environmental changes can be summarized as follows.

4.1.1. Late glacial to Early Holocene deposits

In general, below the lower part of the Holocene sediment lies the motley-colored clay of the late Pleistocene which resulted from weathering conditions at the LGM, the exception being areas where the incised valleys and post-LGM paleo-channels show erosion and downcutting of the pre-existing mottled Pleistocene sediment layer. Based on the thickness of Holocene sediments, isobaths mostly lie below -10m, deepening from the northwest to the southeast and reaching their greatest depth (about -40m) in Lingding Bay.

Early Holocene sediments in the Pearl River Delta can be dated back to about 11000 years ago [38]. In most cores, sediments from the early Holocene (11000-9000 years BP) are mainly composed of coarse sand, indicating a fluvial sedimentary environment. The average particle size from the key boreholes is about 125-250µm, with medium fine-sand being dominant, making up an average of 62% of all particles. Gravelly-sands from paleochannels can be found at the bottom of a few boreholes, while organic-rich layers are more common in most cores. Such organic-rich layers, mainly comprising peat deposits with macro remains of plants, suggest that freshwater wetland may have been more extensive than river channels before this became an estuary environment. This sedimentary stage represents the earliest Holocene deposits in the Pearl River Delta and show gradual changes in the sedimentary environment from fluvial to marine facies.

It should be noted that, with postglacial sea-level rise, some areas in the southern part of the Pearl River Delta with deep incised valleys may already have been inundated by sea transgression around 9000 years ago, possibly even earlier [39]. Sea-level data in the PRD indicates a rise of relative sea level from -49.3 to the present level between 10500 and 7000 cal a BP [14]. One record, from a location in the head of the PRD plain near the confluence of the North and West Rivers, exhibits marine sediment at -23.47m (dated to 9482 years BP) which could be the result of sea transgression along a deep incised valley.

4.1.2. Mid-Holocene marine deposits

During the Holocene marine transgression in the PRD, deeper water sediments of marine facies containing abundant microfossils, such as mud and silt, were deposited over a previous sediment stage which is characterized by shallow marine or fluvial sediments, such as silty or coarse sand. The sediments of the Mid-Holocene (8000-5000 years BP) mainly consist of gray-black clay, silt, and sandy silt, which are often deposited in a tidal flat environment. The abundance of ostracods, foraminifers, saltwater diatoms, and mangrove pollen indicate a marine sedimentary environment. The average particle size of the sediment is about 32-63µm, and fine sand makes up about 50% of the material. This stage represents a period of marine transgression, during which relative sea level in the PRD reached its highest point (or peak transgression) at about 8-7 ka (Figure 3). The results of 14C dating from DH7 show that the interval with the maximum content of foraminifera and the highest Sr/BA ratio is between 7984-7441 years BP, suggesting that this was when the maximum marine transgression occurred. The spatial distribution of sediments associated with the marine transgression extends up to the northwest of the delta, close to the confluence between the West and North Rivers. After 7000 years BP, relative sea level in the PRD changed very little [14]. After the maximum transgression there is more evidence of stronger river input, characterized by fluvial deposits. However, fluvial sediment facies vary spatially due to the complex topography of the paleo-estuary and the different rates at which they were laid down in different places.

4.1.3. Increasing alluvial deposits during the Late Holocene

The lithofacies of the late Holocene, since 5000 years BP, represent a sedimentary stage of marine regression characterized by increasing alluvial sediments which indicate progressive deltaic progradation [11,18,40,41,42]. The characteristic sediments of this stage are composed of silt, fine sand and occasionally coarse sand. The topmost layer is usually yellow-brown clayey silt, indicating the alluvial and fluvial deposits of the deltaic plain setting. In general, the particle size of the late Holocene sediment is relatively fine, about 4µm on average, with silt accounting for more than 50%. The sediments in central and southern parts of the PRD still contain some foraminifera, but the number of fossils is significantly lower. Marine shells are occasionally found in the lower sediment layers of this stage, but some very dense shell deposits may be Neolithic shell mound in relation with human activities. The magnetic susceptibility of this layer generally increases rapidly. The content of herbaceous pollen (especially Poaceae) also increases greatly, together with the dominance of C4 plants inferred by δ¹³C. The average CIA chemical alteration index is about 60-70, which is significantly higher than the 30-40 found in the Mid-Holocene marine sediments. Such characteristics reveal that the sedimentary environment changed from marine to fluvial facies. The accretion sedimentary sequence during the late Holocene reflects increasing fluvial transportation and accelerated river discharges.

4.2. Holocene sediment isobaths and underlying deep incised-valley

In this study, the thickness maps of the late Quaternary and Holocene deposits were drawn based on the core dataset. The geomorphological features of the basement revealed by the Quaternary isobath map are not only fundamental to the sedimentary evolution of the PRD Plain, but also shed light on channel avulsion and the process of marine transgression and regression. According to the spatial distribution of the Quaternary deposits in the Pearl River Delta (Figure 2c), the topography of the delta basement is relatively gentle and the average thickness of the Quaternary deposits is ca. 25 m. In places the depth of the Quaternary deposits reaches -60m, for instance in the Wanqingsha area near the Jiaomen and Hongqimen outlets. There are three bead-like depocenters along the west and north orientation which, together with the bedrock hills, form parallel ridge-valley landscapes. Although these findings are similar to those reported by Huang et al. [7] and Chen et al. [8], it is worth noting that we successfully reconstructed the most accurate Quaternary isobath map so far and obtained a clearer understanding of the spatial distribution of Quaternary deposits in the Pearl River Delta plain by using the borehole database.

The Holocene isobath map indicates that Holocene sediment thickness varies across the region because of underlying bedrock/islands, incised-valley geomorphology, and Pleistocene sediment distribution. Based on geomorphic features and underlying stratigraphy, the PRD can be separated into three zones of differing sediment thickness. The first, with sediment thickness of >25 m, is an elongated zone stretching from northwest to southeast, roughly extending from the west of Foshan to the north of Lingding Bay. This zone reflects the location of the deep incised valleys. The second zone comprises the transitional areas from the deep valley, with a larger surface than the first zone. The third zone, with a thickness of less than 10 m, is located primarily in the northwest of the inner delta area, areas surrounding Panyu terrace, and marginal delta areas along the rivers. Holocene sediments generally thicken towards the southeast, possibly due to the southward dip of the underlying strata.

The isobaths of the Holocene deposit show several zones which are relatively low in altitude and elongated in surface shape. These likely indicate the location of incised valleys with two valleys roughly aligned and growing deeper downstream of the West and North Rivers. The largest incised valley can be tracked upstream to a point of origin northwest of Foshan City. These NW-SE trending fluvially eroded and elongated paleotopographic deep valleys are larger than a single channel, and join together at the lowest seaward reach of the delta to form a broad valley more than 10km wide. The lower end of the larger deep incised valley coincides with the present-day river outlets of Hongqimen and Hengman in the north of Lingding Bay. Another deep river valley originated from the south of present-day Jiangmen city, through an area of numerous islands to the Modaomen outlet and the sea. However, this incised valley, in the lower reach of the current East River, was less significant. The current East River basin only accounts for 5.96% of the entire drainage area of the Pearl River Basin. We suggest that the small water flow of the East River during the last glacial maximum made it unable to erode into a large and deep incised valley.

Formation of deep incised valleys is generally associated with tectonics, river capture, and drops in sea-level [46]. The most accepted and widespread paradigm for explaining the downcutting of incised valleys during the Quaternary is the fall in sea level driven by the increase in global ice volume [47]. In the case of the PRD, this paleo-topographic feature formed as result of a fall in sea level during the LGM and fluvial downcutting on the exposed Pleistocene deposit plain. The valley then began to fill as a result of postglacial sea level rise, accumulating deposits through subsequent Holocene sea-level cycles. Most Holocene sediment profiles in the incised valley in the center of the PRD are characterized by backstepping fluvial and estuarine deposits, overlain by fluvial deposits. In general, there is a surface of unconformity between the areas based on the late Pleistocene layers or bedrocks and the areas of the incised valleys which comprise alluvial sediments followed by estuarine deposits from the time of maximum transgression.

Our reconstructions show that the Holocene beaded depocenters were located in an elongated NW-SE oriented area, rather than in the middle of the delta or at the mouth of the paleo-estuary, a finding which differs slightly from the distribution pattern established by Wei et al. [12]. It has been found that there is a close relationship between tectonic faults and the deep incised valleys. Two NW-SE faults (Beijiang fault and Baini-Shawan fault) control the extent and direction of the valleys. The strip from the West River to the mouths of Hengmen-Hongqimen-Jiaomen is the largest zone of the deep incised valley. Another important incised valley which ended in Modaomen was also formed in association with a long extended fault (Xijiang fault) that shapes the western boundary of the delta. During postglacial sea level rise, the incised valleys in the PRD were all submerged by transgression. In the late Holocene, with the gradual reduction in the estuary’s accommodation space, the delta plain continued to expand, gradually forming the present landform and river network system.

5. Discussion

The spatiotemporal variations in sedimentation rates in the PRD are significant because of the complexity of the estuary topography, spatial heterogeneity of accommodation space, and underlying river valleys. According to the statistical results for sedimentation rates, the mean value for the PRD during the early Holocene is relatively small [7,12]. This can be explained by relatively low sea levels, with most of the delta area remaining exposed prior to 9000 year BP [10,48,49] and a limited marine transgression from south to north which was confined to the deep valleys. Only the low-altitude areas of the southeastern PRD were initially affected by the marine transgression, and here the sedimentation rate accelerated.

Relative sea level rose rapidly from the early to mid-Holocene, from about -49m to the present level between 10,500 and 7000 cal. years BP, which averages at around 13.7mm/a [14]. Sedimentation rates during the early Holocene in the PRD region were relatively small because most of the sediment was transported into the offshore areas in the South China Sea [12]. A substantial volume of proxy data on the marine sedimentary environment indicates that the age of marine sediments may be as high as >9 ka, and the Holocene marine transgression reached its maximum at around 8-7 ka BP (Figure 3), while the coastline of the large paleo-estuary was directly bordered by mountains, hills, and terraces. The broad deltaic basin during this period provided a large accommodation space for the delta plain to subsequently prograde [49]. The sedimentation rate during 7.5-5.0 ka was about 2.34 mm/a on average in the deltaic depocenter [12]. In this study we found that, according to the data from 22 boreholes, and consistent with previous results, the average sedimentation rate for the PRD as a whole was 2.14 mm/a in 7-6 ka, and then gradually decreased between 6-5 ka (1.91 mm/a) and 5-4 ka (0.94 mm/a). The reduced sedimentation rate was considered to be the result of a progressive decrease in monsoon-induced fluvial discharge [6]. However, the sedimentary rate calculation did not take into account the hiatuses (time gaps in the sediment sequences, widespread distributed in PRD) in the Holocene stratigraphic record caused by strong erosion during the late Holocene. In general, erosion tended to occur when there was a lack of accommodation space. The results of dating the material in our dataset show a gap in the Holocene sediment record between about 8-4 ka, illustrating a hiatus in the stratigraphic sequence (Figure 5). This may be the result of sediment erosion by stronger river runoff during the late Holocene. At the beginning of marine regression, coastal sedimentation was largely ascribed to climatic factors [49,50], and partly to regional differential vertical crustal movements [51]. The Mid-late Holocene transition was a period of profound changes in summer monsoon and hydrological conditions in China. A great number of studies have shown that the runoff of the West River and North River was significantly stronger in the late Holocene, since ca. 4000 years BP, coinciding with the decrease in accommodation space. Increasing river runoff and flooding led to the large-scale erosion of underlying early-mid-Holocene sediment layers and/or even Pleistocene deposits, particularly in areas of high hydrodynamic energy. In addition, several studies concluded that the sedimentary characteristics of this period were caused by the unique and complex geomorphic structure of the Pearl River Delta, which was subject to complex and heterogeneous sedimentary processes [12,28].

In East and Southeast Asia, significant delta progradation after the mid-Holocene is considered to have resulted from forced regression [49,52,53]. Transport of sediment from land to the delta via rivers, land erosion, and soil loss increased progressively. Existing research offers two main schools of thought on the Holocene changes in the coastline of the Pearl River Delta. The first is that, since the maximum marine transgression at about 8-7 ka, the coastline of the Pearl River Delta advanced sequentially from the head deltaic areas (topsets) towards the sea in a roughly parallel expansion pattern (Figure 6). This pattern of coastline changes is supported by archaeological findings, historical archives, family ancestral books, and modern survey data spanning from the Neolithic Ages through the agricultural period to the modern industrial era [6,10]. This point of view is common and generally accepted by geologists and geographers. The alternative point of view is that the evolution of the coastline of the PRD was not driven from north to south, or from the river outlet to the downstream area. Rather, multiple deposition bodies around the islands or sand mounds simultaneously developed in the north and south without overlapping. This perspective emerged from simulation results based on hydrodynamic and bed topographic data using the Pearl River Delta Long-term Morphodynamic Model (PRD-LTMM), reported by Wei et al [39,54] and Wu et al. [55]. They suggest that there is a linear relationship between energy flux and cumulative deposition volume, and that energy flux decreased more in the outer part of the estuary than in its inner areas. We believe that both of the above perspectives have a sound basis in theory and evidence. However, a new study has pointed out that preferential sedimentation alone was not sufficient for the emergence of delta plains around the islands [6]. Therefore, on the basis of the reconstruction of isobaths conducted for this study, we favor the coastline change pattern proposed by Zong et al. [10] and Xiong et al. [6], that the seaward delta progradation was roughly parallel to the paleo coastlines. Only a few areas of the deposition bodies around the islands [12], tested through sediment lithofacies, showed signs of adaptation. It has been suggested that the rate of sedimentation close to the islands slowed significantly when the accommodation space was reduced to just a few meters, but only when the mouths of distributaries reached the islands, sedimentation resumed and delta plains emerged [6,56]. Therefore, we speculate that some of the sedimentary bodies around bedrock islands simulated by the hydrodynamic model may be somewhat exaggerated in terms of surface size. However, we adapted the basic framework simulated by the model (PRD-LTMM) [54], which is determined by water energy flux and its spatial variation, to reconstruct the paleo channels.

Although humans arrived in the Neolithic age (some 7-6 ka) [57] neither the climatic factor nor early human activity, including the shift from hunting and gathering to early agriculture around 5.0-4.5 ka [58] had a significant impact on the formation of the delta. Our calculation of the average sedimentation rate in the PRD over the past 4000 years shows a gradual increase, reaching 3.61mm/a between 2-1 ka and 5.24mm/a between 1-0 ka, the highest rate in the Holocene. This is consistent with previous findings that the average sedimentation rate during the development of the PRD has been at its highest since 2500 BP, because frequent human activities have increased the amount of sediment discharge to the delta [6,7,15]. However, according to Wei et al. [12], the sedimentation rate of the overall delta has not necessarily increased during the last 2500 years. With the delta plain prograding seaward, accommodation space reducing, and paleo channels narrowing, a large proportion of sediment discharge is transported directly into the lower part of the PRD. In fact, the calculation of the average sedimentation rate in the delta depends on the location of the collected boreholes. Most of the borehole data with dating information in this study are derived from the middle and lower reaches of the delta.

Comparison of historical population evolution [59] with our reconstruction of the Holocene sedimentation rate, coastline changes, and quantitative land growth in the PRD suggests that human activity has played an important role in the evolution of the delta plain. Since the Qin-Han Dynasties (around 2000 years BP), humans have transformed the transport and deposition of sediment in the deltaic areas where human activity has dominated. On the one hand, human slash-and-burn farming and other agricultural activity since the Qin Dynasty led to forest destruction and soil erosion, resulting in accelerated river sediment discharge and delta plain expansion. On the other hand, the newly increased delta land area created more space for people to inhabit, leading to further population growth in the region (Figure 7). Our data provides a mechanism for quantifying anthropogenic changes to deltaic sedimentary systems and topographies.

The particularly high rate of land growth and rapid coastline shift in the past few hundred years has been interpreted as showing that human activities have significantly increased fluvial sediment delivery. Since the Ming Dynasty, wide and scattered sandy lands gradually accumulated in the area of the southern central Pearl River Delta close to the paleo shorelines, and the river network formed during the late Ming and Qing dynasties. The comparison between modern rivers and deep incised valleys shows that the major watercourses today (especially the lower West River) are close to the paleo-channels of the incised valleys. That is, the evolution from paleo to modern river courses has a certain inheritance (Figure 8). However, the complex network river system formed through the interaction of nature and human activities over the past few hundred years. The distribution of Holocene riverbed sand bodies [60] which our results illustrate show that an environment characterized by scattered sandy lands, which play a role in the formation of the modern river network, emerged between 600 and 300 years ago. That is, the modern watercourses are generally located beside the Holocene sand bodies. This landscape of scattered “sand land” surrounded by water was transformed from the Song Dynasty onwards into lands for rice cultivation, with easy irrigation and, subsequently, fish ponds.

6. Conclusions

For this paper, we established a borehole dataset which contains information from over 2800 sediment cores in the PRD, the highest resolution dataset so far available. This dataset makes it possible to demonstrate the spatiotemporal characteristics of the marine and terrestrial stratigraphy and sedimentary evolution history of the PRD. High-resolution isobath maps of Quaternary deposits were created, providing brand-new insight into Holocene sediment distribution, showing a beaded and elongated zone of depocenters (generally >20 m), which mainly run from the northwest to the south of the delta. This reflects that the Holocene basal topography, with deep incised valleys forming during the LGM, was largely controlled by the faults which run in a NW-SE direction. Two subparallel deep incised valleys increase in both size and depth towards the downstream of the West River and North River. The largest incised valleys may reach 10 kilometers in width near Hongqimen and Hengmen, and can be traced upstream to a point of origin to the northwest of Foshan City. We also found that there was no obvious deep incised valley downstream of the East River possibly due to small water flow during the last glacial maximum.

Multiple proxy data from the cores in the PRD show that the sediment interval with the greatest abundance of marine indicators (ostracods, foraminifers, saltwater diatoms, and mangrove pollen) is concentrated within the period between 8000 and 5000 years BP, and the Holocene marine transgression reached its maximum at around 8000-7000 years BP. Sedimentary facies of the late Holocene since 5000 years BP is characterized by increasing alluvial sediments, revealing a rapid transformation from marine to fluvial facies. Around 2000 years BP, in the context of little accommodation space in the PRD, substantial river discharge and increasing fluvial transportation accelerated the evolution of the deltaic plain.

Based on the dating results from 22 boreholes and the reconstructed paleocoastline, we calculated the average sedimentation and land growth rates since 7000 years BP. The reduced sedimentation rate during 7000-4000 years BP is likely to result from a large sedimentary accommodation space and progressive decrease in the monsoon-induced fluvial discharge. Our calculation of the average sedimentation rate over the last 4000 years shows a gradual increase, with the rates between 2-1 ka and 1-0 ka BP reaching their highest during the entire Holocene. Considering with historical population evolution and quantitative land growth of PRD, we conclude that human impact became increasingly important to the evolution of the delta plain [6,7,12,15]. The comparison of major modern rivers and reconstructed deep incised valleys demonstrates that the evolution from paleo to modern major watercourses (especially the West River) has a certain inheritance. Finally, our results regarding the distribution of the Holocene riverbed sand layer [60] together with the landscape of scattered “sand land” surrounded by water which has characterized the area over the past few hundred years since the Song Dynasty, explains the complexity of the modern network river system as a result of the interaction between natural forces and human activities. As well as new data, this study provides a better understanding of delta evolution that can be extended to the assessment of future coastline scenarios and may be used as a basis for coastal zone management and urban planning.