Sampling surface particle size distributions and stability analysis of deep channel in the Pearl River Estuary

Particle size distributions (PSDs) of bottom sediments in a coastal zone are generally multimodal due to the complexity of the dynamic environment. In this paper, bottom sediments along the deep channel of the Pearl River Estuary (PRE) are used to understand the multimodal PSDs′ characteristics and the corresponding depositional environment. The results of curve-fitting analysis indicate that the near-bottom sediments in the deep channel generally have a bimodal distribution with a fine component and a relatively coarse component. The particle size distribution of bimodal sediment samples can be expressed as the sum of two lognormal functions and the parameters for each component can be determined. At each station of the PRE, the fine component makes up less volume of the sediments and is relatively poorly sorted. The relatively coarse component, which is the major component of the sediments, is even more poorly sorted. The interrelations between the dynamics and particle size of the bottom sediment in the deep channel of the PRE have also been investigated by the field measurement and simulated data. The critical shear velocity and the shear velocity are calculated to study the stability of the deep channel. The results indicate that the critical shear velocity has a similar distribution over large part of the deep channel due to the similar particle size distribution of sediments. Based on a comparison between the critical shear velocities derived from sedimentary parameters and the shear velocities obtained by tidal currents, it is likely that the depositional area is mainly distributed in the northern part of the channel, while the southern part of the deep channel has to face higher erosion risk.


Introduction
Particle sizes and their distributions are the fundamental properties of sediments. They have a profound influence on many other properties, including suspension, flocculation, settling, deposition, erosion and re-suspension (Blott and Pye, 2012). In general, sediment particle sizes display a unimodal distribution controlled by a single factor in a certain transport medium. In a coastal region, suspended particulate matter is composed of various constituents including clay minerals, organic matter, salts, trace metals and so on (Maggi, 2009). The mixture of clays/silts and sands will continuously change their relative mass fractions and alternately develop a unimodal, bimodal or multimodal particle size distribution in the flow varying with the tidal cycle in the coastal zone (Yuan et al., 2009). As the multimodal particle size distribution is a routine procedure in many sci-entific and engineering investigations to understand sediment sources, migration trends and marine dynamics (Blott and Pye, 2012), the environmental interpretation of sediment particle size distribution is a frequently discussed topic in sedimentology (Sun et al., 2004).
Usually, hydrometer/sieve, scanning electron microscopy (SEM), X-ray attenuation and laser diffraction techniques are used to quantify particle size distribution (Cheetham et al., 2008). Since sediments are usually multimodal, the division of individual components in bulk samples is necessary. However, it is found that sediments often consist of overlapping components that should be partitioned for the investigation of sedimentary dynamics. A graphical method has been introduced in recent years to distinguish the overlapping components and relate a particular cumulative log probability curve shape in a specific envir-onment. Several processes are involved in the formation of sediment and each process is assumed to generate subpopulations, which separate as lognormal particle size distributions. In mathematical terms, a multimodal sediment can generally be described by the sum of some lognormal functions (Påsse, 1997). The decomposition of observed sizedistributions has also been used to characterize sedimentary facies and depositional environments (Bartholdy, 2006). As a result, they are widely applied to analyze a variety of particle matters in different regions (Bartholdy et al., 2007;Fettweis et al., 2012;Lee et al., 2012).
Previous researches in this area often focus on the source, path and flux of suspended sediment (Zhang et al., 2013), erosion and siltation (Heise et al., 2010), as well as flocculation and settling features of suspended sediment (Xia et al., 2004). These studies take the multimodal particle size distribution as a single concept and use the statistical parameters to describe the whole distribution, such as mean, sorting, and skewness. In later studies, the spatial patterns of summary statistics were assumed to reflect the net sediment transport pathways (Gao et al., 1994;Le Roux and Rojas, 2007). However, few researches have been done on the multimodality of particle size distributions of the bottom sediment in this region. The particle size distributions (PSDs) are regarded as unimodal, while the PSDs are often bimodal or multimodal. Accordingly, it is of great significance to decompose multimodal PSDs into lognormal subpopulations. Through this method, the better understanding of the sediment distribution as well as the complex sedimentary dynamic environment could be obtained. Therefore, the objectives of this paper are, (1) to use a statistical method to divide observed PSDs into log-normal subpopulations, and (2) to investigate the interrelations between the tidal dynamics and the particle size of sediment for better understanding of the potential erosion in the deep channel of the Pearl River Estuary (PRE).
The PRE (Fig. 1), located in a transitional zone between the Pearl River and the South China Sea (lat. 21°20'N to 23°30'N and long. 112°40'E to 114°50'E), is a bell-shaped estuary, approximately 4 km wide at the upstream end (Humen) of the bell and nearly 60 km at the downstream end (between Lantau Island and Macau). The estuary is topographically characterized by three shoals (water depth<5 m) and two channels (water depth>5 m). From the east to the west, the shoals and channels are the East Shoal, East Channel, Middle Shoal, West Channel and West Shoal, respectively. The West Channel, which is also called the Lingdingyang Channel, runs down the PRE and connects with the Lantau Channel (Xia et al., 2013). As the Lingdingyang Channel has been dredged continually, we mainly focus on the PSD in this channel in this paper.

Methods
The 2-cm thick surface layer of the bottom sediment samples were collected at 20 sites (S1-S20, Fig. 1) from the coast to the South China Sea basin in August 2007. Samples were taken from the top section of each grab sample to analyze the grain-size distributions. Grain size analyses of the bottom sediment samples were carried out by using a CILAS 940L laser grain size analyzer (for sediments finer than 1.4 mm) or by sieving (for sediments coarser than 1.4 mm). The Acoustic Doppler Current Profiler (ADCP) was also employed to measure the tidal currents. Meanwhile, a measuring network was setup along the Lingdingyang Channel, which has lasted for about 30 hours in the same time with the sediment sampling (Stations A-D, Fig. 1). The observed PSDs are divided into the log-normal subpopulations by using the commercial fitting program DistFit TM (DistFit TM , 2004). The best quality PSD is defined by the minimum error between the simulated and measured PSDs (Whitey, 1978).
where V is the volumetric concentration; D is the particle diameter; dV and dD are the partial derivatives of V and D, respectively; σ i is the geometrical standard deviation; is the geometric mean diameter; and is the volumetric fraction of the i-th unimodal PSD. The volumetric concentration of an aggregate size interval ( ) is obtained by Eq.
(2), and σ i can be directly obtained from the size distribution plotted by DistFit TM . Six parameters for , σ i and are required to fit Eq. (1), dV/dD is the normalized volumetric fraction by the width of the size interval and is used for the curve fitting to a lognormal distribution (Hinds, 1999). For two modal peaks, the standard deviation of a subordinate PSD (σ i ) is allowed to vary between 1 and 2.5 to prevent unrealistically wide PSDs. The choices of the parameters are based on assumptions and experiences .

Sediment components and the PSD features of the PRE
Based on the classification of Shepard (1954), the median grain size μ(mm), sorting coefficient σ, and skewness Sk, were calculated by using the statistic moment method (Mc-Manus and Duck, 1988) where i is the size class, P i is the percentage of size S i , and n is the total number of size classes. The median size is calculated from samples. The bottom sediments along the Lingdingyang Channel can be classified as silty clay, clayey silt, and sand-silt-clay. Of which, silty clay and clayey silt are the dominant sediment types, 55% and 35%, respectively (Table 1). In addition, sand-silt-clay (10%) is found between the Neilingding Island and the Lantau Island. Table 1 shows the statistical parameters of the bottom sediments along the longitudinal of the Lingdingyang Channel. It is obvious that the median sizes on the surface of the PRE range from 0.003 mm to 0.0344 mm. The median sizes in the northern estuary are relatively fine, with the values ranging from 0.0033 mm to 0.0098 mm. To the south of the Neilingding Island, the bottom sediment is also fine with the median size of about 0.012 mm. The sorting coefficients of the bottom sediments in most areas are larger than 2, while in the northern estuary, the minimum value is 0.59. Furthermore, the sorting coefficient reaches its maximum value of 3 in the middle of the Neilingding Island and near the Lantau Island (S11). These relatively high values indicate that the sediments in this region are very poorly sorted. A clear spatial trend of the sorting coefficient is not found. Positively and negatively skewed sediments are found in the Lingdingyang Channel, with values varying from -0.85 to 1.08. The skewness seems highly correlated with the median size (R 2 =0.75), which indicates that the skewness changes from positive to negative when the particle size decreases. Particularly, the median size, 0.006 mm, is roughly the division between the positive and negative skewnesses. Particle size distributions reveal the positive skewness with relatively high values between 0.03 and 1.08, while the negative skewness with relatively low values between -0.03 and -0.85.

Decomposition of the PSDs and modeling subpopulations
A mathematical method is used to analyze the compositions of the bottom sediment samples in the Lingdingyang Channel. Almost all samples yield acceptable partitioning results with a relatively small residual (2.86%-9.63%). The observed multimodal (2-peaked) PSD is assumed to be formed by overlapping two unimodal PSDs, which implies that the bottom sediments in the PRE have two basic sources. The surface sediment is mainly originated from river runoff and tidal currents in the PRE (Jiang and Zheng, 2008). In addition, the median sizes of the former two sub- populations are very similar when the multimodal PSDs is decomposed into three subpopulations (not shown). The function used in this paper is a sum of two lognormal functions, which is suggested to be a more appropriate function for defining the PSDs. The lognormal distribution could rationally describe the characteristics of the distribution and the causes of the formation of each component of the sample. Its parameters also have definite physical meanings that can reasonably perform the PSD characteristics of each component. Fig. 2 shows the partitioning results of several sediment samples at each station. The bottom sediment samples in the deep channel are composed of fine components with the median size range of 0.0056-0.0065 mm, overlapping with relatively coarse components. It can be seen that the standard deviations range from 1.623 to 2.234, and show little spatial variation ( Table 2). The average standard deviation of the fine components is 1.84, indicating a little poor sorting in the fine components over the deep channel.

Particle size distributions
The relatively coarse components have the median diameter range of 0.0683-0.0928 mm. The average median diameter of the coarse components is 0.0812 mm with little spatial variation. The standard deviations range from 1.682 to 1.985 and the average value is 1.907, which indicates a poorer sorting than the fine components. The change in volumetric fraction of the relatively coarse components is opposite to that of the fine components. The relatively coarse component consists of the main part of the sediment from the upper part to the lower part of the PRE.
The spatial variations in fine and relatively coarse components show that the bottom sediment samples in this deep channel are composed of fine components with the median size range of 5.6-6.5 μm, overlapping with relatively coarse components.
The distribution of the averaged particle sizes for the 20 stations was divided into two lognormal distribution sub- populations, which represent the fine and relatively coarse components (Table 2 and Fig. 3). The averaged median size of the fine components is 0.0058 mm and the standard deviation is 1.771, indicating a little poor sorting in the PRE region. For the relatively coarse components, the averaged median size is 0.0812 mm and the standard deviation is 1.895, which suggests poorer sorting than that of the fine components. Overall, due to the relatively small value of volumetric fraction of the fine components, the change in the volume of each component is not obvious from north to south in the deep channel of the PRE.

Discussions
In this section, the sediment transport and hydrodynamic parameters are used to evaluate the possibility of erosion and siltation in the Lingdingyang Channel.

Sedimentary dynamic processes
The PRE is mainly characterized by the combined ef- fects of the open-sea tidal currents and large runoff from the Pearl River (Ji et al., 2011a(Ji et al., , 2011b. To the north of the Neilingding Island, the currents are rectilinear, while in the south of the island, the currents are rectilinear in the channel and rotary in the shallow waters. The result of tidal current measurements shows that the average flood and ebb tidal current velocities are 0.465 m/s and 0.690 m/s, respectively at Station A, which implies that the ebb current is stronger than the flood current. However, the maximum velocity of ebb current at Station A is smaller than 1 m/s (Table 3), which may make the fine components in the suspended sediment deposit easily. The result also indicates that relatively larger amount of fine components in the bottom sediment are around this area. Station B is located in the maximum turbidity region of the PRE, which is characterized by high suspended sediment concentration. The field measurement shows that the maximum ebb current is over 1 m/s (Table 3), which may enhance the sediment transport ability and wash away fine components. The maximum velocities of the flood and ebb currents are both close to 1 m/s. However, the average flood velocity is larger than the average ebb velocity. This difference may also increase the sediment transport ability and wash fine components landward. It can be noticed that the velocities in flood and ebb are both smaller than 0.6 m/s at Station D, located out of the estuary, which probably makes the suspended sediment tend to deposit. The results of a numerical flow model in the PRE are employed in this paper to explain the tidal current variations in this area. The model system using a modified version of the Princeton Ocean Model (POM), which is driven  by tides, meteorological forcing and buoyancy forcing, associated with freshwater runoff from the Pearl River system (Ji et al., 2011a). We focus on the 3D estuarine circulation over the PRE and adjacent coastal waters simulated by the fine-resolution inner model (Fig. 4). As the currents in the wet season were calculated by this model system over a 3year period, the flow patterns represent well the flow behavior in the wet season.
In the northern and central parts of the Lingding Bay, the mean currents at both the surface and bottom were directed southward toward the mouth of the estuary, indicating that the river discharge created a dominant flow for the nontidal current in the estuary (Fig. 4). The bottom currents in the outer Lingding Bay and adjacent coastal waters are mainly landward. Larger bottom currents appear in the offshore deep waters (Fig. 4b) (Ji et al., 2011a). The main flows show that the current at the bottom was dominated by the seawater intrusion while the surface is dominated by the river discharge into the Lingding Bay.

Potential erosion, siltation and channel stability
Most of the sediment erosion or transport formulas are developed for uni-directional flows. These formulas application is restricted to tide-dominated oscillatory flow and tidal currents are considered to be quasi-steady (Hayter et al., 1990). In order to analyze the scour and accretion on the seabed in Lingdingyang Channel, the shear velocity in the boundary layer and the risk of erosion are evaluated in this paper. Moreover, the Karman-Prandtl equation is used to calculate the shear velocity.
In Eq. (6), U 0 and U(z) are the shear velocity (m/s) and the current velocity at z above the ground (m/s); z is the height of a current point above the ground (m); z 0 is the bed roughness length (m) and it reflects the total roughness length, which is used to calculate the total bed shear stress. The value of z 0 ranges from 0.07 cm to 0.25 cm (Drake et al., 1992) while the global method gives a mean z 0 value of 0.134 cm (Cheng et al., 1999). In this paper, the main value for z 0 is assumed to be 0.16 cm. κ is the von Karman constant (0.4).
In this paper, the critical shear velocity is used instead of the critical shear stress. According to Zanke (1977), Eq. (7) has been extended by a term which describes the cohesive and viscous influence for finer material.
where u cr is the critical current velocity (m/s), ρ′ is the relative density, g is the gravitational acceleration (9.81 m/s 2 ), d is the grain size median (m), and ν is the water kinematic viscosity (m 2 /s). In this paper, ν is assumed to be 0.8962× 10 -6 m 2 /s according to the underlying average temperature of 25°C in summer. And c is the coefficient for the critical velocity depending on the level of consolidation. Based on the "quadratic friction law" (Poulos and Chronis, 2001), the critical shear velocity can be calculated from the critical current velocity in Eq. (8): where is the critical shear velocity (m/s), C D is the drag coefficient, and C D (z, z 0 )=[κ/ln (z/z 0 )] 2 . By combining Eq. (7) with Eq. (8), the critical shear velocity can be calculated by where ρ s and ρ w are the density of grains (2650 kg/m 3 ) and the water density (kg/m 3 ). The risk of erosion (R) will be calculated by the empirical Eq. (10) introduced by Zanke (1990).
where τ * is the bed shear stress (N/m 2 ); τ cr is the critical bed shear stress (N/m 2 ). The value of R ranges from 0 to 1, the  erosion risk probabilities become higher when the value of R increases. The patterns of the erosion risk show a higher percentage for areas with the median of fine sand in the Lingdingyang Channel. The value of R is calculated from the simulated and measured data in this paper. Consequently, if erosion occurs, the risk of erosion (R), which is dominated by the grain size or morphological settings, often shows a high value (erosion risk≥90%). Fig. 5 shows the relationship between the tidal currents (shear velocity) and the bottom sediment (critical shear ve-locity) in the measured and simulated data during spring tide in the wet season. The shear velocities that calculated from the measured and simulated data are similar. The critical shear velocity is similar in most parts of the deep channel because of the similar PSD of sediments. Station D is not used to analyze the potential erosion of the Lingdingyang Channel because it is located in the open sea area. In addition, Station B is also excluded from the stability analysis as the high suspended sediment concentration may have great influence on the accuracy of the research result.
It can be seen clearly from Fig. 5a that the critical shear velocity at Station A, which is located at the northern part of this channel, is around 0.04 m/s. This value is always larger than the shear velocity during the measuring time in the wet season. And the erosion trend of the simulated data is similar to that of the measured data. This implies that the erosion probably does not occur around this area. However, it is supposed that sediment gets deposited at the time of slack water because the sediment re-suspension and deposition are highly related to the tidal velocity. For Station E, located at the central part of the Lingdingyang Channel (Fig. 1), the shear velocity is smaller than the critical shear velocity in most of the measured time and this situation is just reversed at the other 15% of the time (Fig. 5b). As for the simulated data, the shear velocity is below the critical shear velocity all the time. The periods that R over 0.9 only occupy a lower proportion of the whole time, which implies that the chance of erosion around this area is small. Station C is located at the lower part of the deep channel (Fig. 1). Interestingly, the shear velocity here is above the critical shear velocity most of the time (more than 70%, Fig. 5c), which means that the erosion probably occurs frequently. Meanwhile, as the period that the value of R is larger than 0.9 makes up more than 50% of the time, Station C suffers a high erosion risk.
According to Fig. 6, the shear velocities are smaller than the critical shear velocities and the risk of erosion is close to zero during the period of the neap tide in the wet season at  FENG Hao-chuan et al. China Ocean Eng., 2017, Vol. 31, No. 3, P. 299-307 Station A as well as Station E. While in Station C, the period that the value of R is larger than 0.9 accounts for more than 40% of the time in neap tide. Fig. 7 presents the comparison of the shear velocity, the risk of erosion and critical shear velocity calculated from the simulated data during the dry season. Apparently, both in Stations A and E, the shear velocities are always smaller than the critical shear velocity in the dry season. While in Station C, the period when the shear velocity is higher than the critical shear velocity accounts for nearly 70% of the time during the neap tide and nearly 60% in the spring tide in the dry season. In addition, more than 55% of the time the erosion risk is larger than 0.9 during spring tide while only 30% of the time that happens in the neap tide (Fig. 7).
Based on the above analysis, there is a high risk of erosion in Station C (more than 90%) both during the dry season and wet season, while the risks in Stations A and E are low (most of the time nearly 0) (Figs. 5, 6 and 7). Overall, it can be predicted that erosion happens in the southern part of the Lingdingyang Channel, while the deposition mainly occurs in the northern part.

Conclusion
Results from the curve fitting analysis of the multimodal PSDs suggest that the near-bottom sediment in the Pearl River Estuary generally has a bimodal particle size distribution and the parameters for each modal can be determined. The PSD of bimodal sediment samples can be expressed as the sum of two lognormal functions, representing the relatively coarse and fine grain-size components. The fine components have a median size range of 0.0056-0.0065 mm. The standard deviations of the fine components in PRE are generally smaller than 2. The median sizes of the relatively coarse components range from 0.0683 to 0.0928 mm, and show little spatial variation in the PRE. The standard devi-ations range from 1.682 to 1.985 and generally do not change obviously. The main material is the relatively coarse component at each station of the Pearl River Estuary.
The potential sediment erosion is investigated through the analysis of the interrelations between the tidal currents and particle sizes of the bottom sediments. The critical shear velocities are similar at most parts of this channel due to the similar PSDs. With comparison of the critical shear velocity with the shear velocity, it can be predicted that in the wet season and dry season there would be no erosion in the northern part of this deep channel. A small amount of erosion is likely to occur at the central part of the channel, and erosion may occur very frequently in the southern part of the deep channel.