Numerical simulation for the effects of waves and grain size on deltaic processes and morphologies

1 Introduction

As a low, flat, and fertile land close to the river and the sea, a delta is rich in biological resources and provides favorable conditions for the development of agriculture, fishery, and aquaculture. Therefore, a delta in a climate-friendly area is usually a densely populated and economically developed area in a locale [1]. While deltas account for only 0.56% of the world’s total surface area, 4.1% of the world’s population live in deltas. The population living in deltas is estimated to increase from 237 million in 2000 to 322 million in 2020. The population living in some deltas, e.g., the Yangtze Delta (Shanghai, China), Zhujiang Delta (Guangzhou, China), Mississippi Delta (New Orleans, USA), and Niger Delta (Lagos, Nigeria), grows by 1.59% annually, which is higher than the world population growth rate of 1.11% [2]. In the case of the Nile Delta, there are 50 million people living in it, and the delta contains large cities such as Alexandria and Cairo.

Ancient delta deposits are highly valued because they possess abundant coal, oil, and natural gas. Because of petroleum geological exploration, many oil and gas fields have been found to be related to the sedimentation of deltas, most of which are large or extra-large. For instance, most of the oil and gas reservoirs in the Burgan oil field in Kuwait, the Maracaibo Basin in Venezuela, and the Bolivar Offshore Field are related to deltaic sedimentation [3,4,5,6,7,8]. Hence, studying the formation and evolution of deltas is of great importance in production and human development.

Great emphasis has been placed on studies of deltas, and research results have been fruitful. Delta morphology is mainly dominated by rivers, waves, and tides [9]. Studies have shown that the grain size of sediment also affects delta morphologies [10,11]. Since 2000, as a result of the improvement of computing power and the development of the fluid dynamics theory, geomorphology simulation using computers has been carried out extensively in geological studies. It has proven to be a feasible method for exploring the formation of river-dominated deltas using geomorphological numerical models based on physical models [12]. According to numerical simulations, shallow basins are more likely to form deltas whose branch channels are formed by sedimentation at mouth bars, such as the Wax Lake Delta in Louisiana, whereas deep basins are more likely to form deltas of the classic bird-foot shape [13]. Simulation results also indicated that sand-dominated deltas are more fan-shaped, whereas mud-dominated deltas are more bird-foot-shaped in planform [14]. According to a three-dimensional hydrodynamic numerical simulation, the hydrodynamic conditions in the upper reaches of the river control the planform of the delta, and the conditions in the lower reaches control the migration of the mouth bar and branch channels [15]. The number of branch channels in a delta is closely related to the flow of the river [12,16]. The hydrodynamic-based delta numerical simulation method has been widely applied in geological studies on deltas and has generated ample research results [17,18,19].

As a result of the complexity of sedimentary simulation and the duration of time required for computation, previous numerical simulations often excluded wave action. Although the sedimentary simulation of river-dominated deltas is already well-developed, there is much work to be done on the sedimentary numerical simulation of wave-dominated deltas.

The Mossy Delta, shown in Figure 1a, is located in Saskatchewan, Canada, with median grain-size D₅₀ = 0.125 mm and P_m : P_r = 0.15 (the ocean energy P_m is equal to the sum of the square of the mean monthly maximum significant wave height and the square of the tidal height difference, the river energy Pr is equal to the product of average water discharge and delta plain gradient scaled by 11) [20]. The Niger Delta, shown in Figure 1b, is located in Nigeria, P_m:P_r = 1.2 and D₅₀ = 0.15 mm; the data of P_m, P_r, and D₅₀ for the two deltas were extracted from Table 1 in Syvitski and Saito [20]. The Pm:Pr values of the two deltas in Figure 1 indicate that the Niger Delta was much more affected by waves than the Mossy Delta. Although they had similar sediment D50 values, their morphologies differed greatly, with one having a semicircular shape and multiple distributary channels (Figure 1a) and the other having an arcuate shape and only a few distributary and straight channels (Figure 1b). If the wave was not considered in the numerical simulation, the numerical simulation in previous studies could perfectly restore the shape of the Mossy Delta [21], but not the shape of the Niger Delta.

Figure 1

Different morphologies of deltas with similar D₅₀ values under the effects of different wave energies. (a) Delta front with finger-like mouth bar. (b) Delta front with a smooth arcuate shape. Deltas were located in Saskatchewan, Canada (54°04′N, 102°21′W) (Figure 1a), and Nigeria (4°27N, 6°6′E) (Figure 1b). Images were obtained from Google Earth.

Figure 2 shows images (from Google Earth) of several deltas affected by similar wave energy. According to the collected data, the deltas in Figure 2 are as follows: (a) median grain size D50≥0.5mm, mean significant wave height 1.2 m. (b)D50=0.2mm, maximum monthly significant wave height 1.5 m. (c) D50 unknown, mean significant wave height 1 m. The data were obtained from Syvitski and Saito [20] and Orton and Reading [11]. Although the deltas were affected by similar waves, with an average significant wave height of about 1 m, their morphologies were not similar. This suggested that in addition to the influence of waves, tides, and rivers, the deltaic morphology may also be affected by sediment grain size.

Figure 2

Different deltaic morphologies under similar wave conditions. Deltas were located at the Bella Coola Valley in British Columbia, Canada (52°22′N, 126°45′W) (a), the Iberian Peninsula, Spain (40°42′N, 0°42′E) (b), and Brazil (10°24′S, 36°29′W) (c). Images were obtained from Google Earth.

To investigate the combined effects of waves and grain size on morphologies of deltas, a new numerical simulation model was developed in this study. Two experiments were performed in this model: (1) with all the other conditions set the same, morphologies of deltas with and without wave action were compared. (2) Under fixed wave energy conditions, the effects of different grain sizes of the sediment on delta morphology were studied. This model quantitatively analyzes the development process of deltas and provides a new method for revealing the distribution of sand sediments.

2 Materials and methods

To fully reflect the process of deltaic sedimentation and establish its sedimentary model, the front-based computational fluid dynamics method was used in this study instead of the traditional behavior-based model, resulting in a so-called process-based model. This method comprehensively considers the influence of hydrodynamics on sedimentary grains, realizes the transport and erosion of sediments, and finally changes the landform under hydrodynamic conditions. The accuracy of this method has been verified by abundant numerical simulations [14,15,22,23,24].

Deflt3D is a hydrodynamic numerical simulation software that was developed by Delft Hydraulics in Delft, the Netherlands. Following the Navier–Stokes equation and adhering to sediment-transport rules and the Law of Conservation of Mass, the software describes the motion of sedimentary grains under hydrodynamic forces by solving the equations concerning the transport, sedimentation, and erosion of sediments. In this study, the simulation calculation process was completed using the Delft3D software, the post-processing of data was realized through C# programming, and the data were three-dimensionally displayed using the OpenGL method.

The transport of cohesive and non-cohesive suspended sediment is calculated by solving the depth-averaged version of the three-dimensionally advection-diffusion equation:

where ci is the mass concentration of the sediment fraction i (kg/m³); u, v, and w are the flow velocity components (m/s); εs,xi, εs,yi, and εs,zi are the eddy diffusivities of the sediment fraction i(m²/s); and wsi is the sediment settling velocity of sediment fraction i (m/s).

The sedimentation velocity of the non-cohesive sediments of fraction i(wsi) was calculated using the grain diameters defined by the user [25], such that

where g is the acceleration due to gravity (9.8 m/s²), sⁱ is the relative density of the sediment fraction i, Dsi is the representative diameter of sediment fraction i, and v_w is the kinematic viscosity coefficient of water. For non-cohesive and cohesive sediments, the erosion and sedimentation of suspended sediments were calculated individually. The exchange of non-cohesive suspended sediments with the materials of the outermost stratum can be seen as a combined function of the erosion of sediments and the sedimentation of suspended matters. Thus, the erosion is expressed as a source term and the deposition as a sink term [25]:

where α2i is a correction factor for sediment concentration, εsi is the sediment diffusion coefficient evaluated at the bottom of the kmx of the sediment fraction (i), cai is the reference concentration of sediment fraction (i), ckmxi is the average concentration of the kmx cell of the sediment fraction (i), and Δz is the difference in elevation between the center of the kmx and Van Rijn’s reference height.

The erosion and sedimentation of cohesive sediments are calculated according to the method proposed by Partheniades [26]:

where S(τcw,τcr,ei) is the erosion step function for cohesive sediment,S(τcw,τcr,di) is the deposition step function for cohesive sediment, τcw is the maximum bed shear stress due to current and waves (as calculated by the wave–current interaction model selected by the user), τcr,ei is the user-defined critical erosion shear stress (N/m2), and τcr,di is the user-defined critical deposition shear stress (N/m2).

The waves were calculated using the Delft3D-wave, which uses the third-generation Simulating Waves Nearshore (SWAN) wave model in this study [27,28]. The wave model was developed based on the second-generation hindcasting of waves in shallow-water wave model [29]. Because of the completely implicit scheme that has been implemented, the wave computation in the SWAN wave model can be unconditionally stable, which greatly improves the computational efficiency. In the SWAN wave model, the evolution of the spectrum is represented by the spectral action balance equation, which is expressed in the Cartesian coordinates system as follows [30]:

where the term N is a function of η and θ and the term N denotes the action density spectrum. In equation (7), η represents the relative frequency and θ denotes the wave direction. The first term on the left side of equation (7) stands for the local rate of change of action density in time and the second and third single terms represent a process of propagating in the x direction at a cx speed in space and a process of propagating in the y direction at a cy speed in space. The fourth term denotes the relative frequency shift caused by changes in water depth and water flow (with propagation velocity cη in η-space). The fifth term represents depth-induced and current-induced refraction (with propagation velocity cθ in θ-space) [31,32,33]. On the right side of equation (9), the term S is a function of η and θ, S is the source term in terms of energy density representing the effects of generation, dissipation, and nonlinear wave–wave interactions. Waves computed in the SWAN wave model, operating as input conditions, participated in computing the migration, erosion, and sedimentation of sediments, ultimately affecting delta sedimentation.

The transport of bedload is calculated by the method described by Van Rijn [34]:

where veff=vR2+Uon2, Sbi is the bedload sediment discharge per unit width of the ith sediment fraction (m2/s), ρs is the specific density of sediment (kg/m3), wsi are non-cohesive sediment settling velocities of the ith sediment fraction, Dsi is the representative diameter of sediment fraction i, vR is an equivalent depth-averaged velocity computed from the velocity in the bottom computational layer, assuming a logarithmic velocity profile (m/s), Uon is a near-bed peak orbital velocity (m/s) in the onshore direction based on the significant wave height, vcr is the critical depth-averaged velocity (m/s) for initiation of motion based on a parameterization of the Shields curve, and si is the relative density of the sediment fraction i.

As tide and deltaic deposits accumulate in the fluvial-to-marine transition zone, the seaward part of estuaries and deltas is subjected to tidal action that produces an alternation of landward-directed and seaward-directed tidal currents. And this area is one of the most complicated areas on earth, which results in a large number of terrestrial and marine processes [35]. Although we can achieve tidal effects by changing the water level and estuary flow direction, more research is needed on the accuracy of the simulation results. So the effects of tide are not considered in this study.

3 Numerical simulation of sedimentation

The same grid was used for all the simulation experiments. The grid consisted of 200 × 150 cells, each of which had a length and width of 25 m, and the simulated computational domain was a 5,000 m × 3,750 m basin with a bed slope of 0.000376 to the south (Figure 3). The thickness of the initial stratum was 5 m, and the grain size of the initial stratum was distributed as the same as the grain size of the sediment input in this simulation. Moreover, the spatially constant Chézy value for hydrodynamic roughness was 45. The water discharged into the basin was Q=1,000m3/s. The cohesive and non-cohesive sediments input into the model totaled 0.1 kg/m³ and remained constant in terms of time. To ensure that the delta can enter a stable state, the simulated duration was about 8 years, and each model was run iteratively 2,10,00,000 times at a time step of 0.2 min. Because sedimentation is slow in a real sedimentary environment, the morphological bed updating factor was set as 175, which not only reduced the number of iterations to 1,20,000 but also maintained the accuracy of the geomorphological evolution [15,23,24]. Moreover, the initial river channel was 250 m wide, 2.5 m deep, and 500 m long. The water surface was 0 m high. The wave action was not introduced in a small number of simulation experiments. As for the other simulation experiments, the waves were spread from south to north; that is, the waves were perpendicular to the initial coastline. In relatively shallow areas (coastal seas), wave action becomes important. So in our numerical experiments, the effects of waves on current and sediment transport (via forcing, enhanced turbulence, and enhanced bed shear stress) are accounted for; they are: (i) forcing by radiation stress gradients [31,36]; (ii) Stokes drift and mass flux [37]; (iii) streaming [38]; (iv) wave-induced turbulence [39]; and (v) enhancement of the bed shear stress by waves. In our study, the basic parameters of the wave [40], based on the actual significant wave heights of the three different deltas in Figure 2, the significant wave height was set to 1 m, the peak period was set as 5 s, the wave angle was 0° (perpendicular to the initial coastline), and the wave parameters remained constant over time. Table 1 shows the basic parameters for all numerical simulations.

Figure 3

Numerical simulation of the initial model. The water surface was 0 m high, the initial river channel was in the north of the model, and the water and sediments entered the system from the north channel. If there were waves, the waves spread from south to north. The red line denotes an open boundary.

For a more realistic simulation of the input sediment, the grain size of the sediment input into the simulation system was assumed to be normally distributed. Since only discrete data can be input into the computer simulation system, the incoming sediment discharge was partitioned into five different sediment fractions during each model run (Figure 4a). By changing the standard deviation (σ in Figure 4b) and median grain size (ϕ50 in Figure 4c), a variety of different distributions of the grain size were created. The grain size (D50) that followed the actual global wave-dominated delta data [11,20] is mainly distributed between 0.15 and 0.5 mm. According to the grain size data and preliminary numerical simulation results, D50 is limited between 0.125 and 0.5 mm. Cohesive sediment was defined as sediment with a grain size ≤64 µm, which is about 3.97ϕ. Hence, changes in standard deviation and median grain size also changed the content of cohesive matter in the sediment. In this study, we used Excel’s random number generating formula for simulating the grain size which is consistent with the standard deviation (σ) formula (Folk, 1974) [41], and these data are used to calculate the properties of sediment grain size (% cohesive, % non-cohesive,ϕ25, and so on). These random numbers have been divided into five fixed grain size to represent the grain size distribution of the input sediment. No doubt, using more sediment fractions creates a more highly resolved grain-size distribution, but sensitivity tests show that a grain-size distribution discretized by five sediment fractions results in a sufficiently accurate result and a tolerable calculation time.

Figure 4

Example grain-size distributions. (a) The continuous grain-size distribution was discretized into five sediment fractions as denoted by the black rectangle, where σ=1ϕ and D₅₀ = 0.25 mm. (b) The ranges in standard deviations, D₅₀ = 0.25 mm. (c) The ranges in median grain sizes, σ=1ϕ.

The Run ID in Simulation Experiment 1 was represented by the capital letter A. To reflect the difference between delta morphology with and without wave action, four experiments were designed for Simulation Experiment 1. The specific parameters are shown in Table 2, where A1 and A3 were not wave-dominated, and A2 and A4 were under the action of waves, and the parameters were the same; that is, the significant wave height was 1 m, the peak period (Tp) was 5 s, and the waves spread from south to north, perpendicularly to the initial coastline (Figure 3). The aforementioned parameters remained constant during the simulation and were mainly based on the actual delta data in Figure 1. The purpose of the simulation was to reproduce the morphologies of the two deltas in Figure 1.

The Run ID in Simulation Experiment 2 was denoted by the capital letter B. In order to reflect the influence of change in grain size of the sediment on delta morphology, nine experiments were designed in this study, the specific parameters of which are shown in Run ID B1–B9 in Table 3. Wave conditions were introduced in Run ID B1–B9 with the same parameters: the significant wave height was 1 m, the peak period (Tp) was 5 s, and the waves spread from south to north, perpendicularly to the initial coastline. In simulations B1–B9, all the conditions were identical except the grain size of the sediment. The main simulation parameters used in these experiments were based on the actual delta parameters in Figure 2. The purpose of the numerical experiments was to reproduce the morphologies of the three deltas in Figure 2.

4 Results

The finite difference calculation method based on the principle of fluid dynamics was used in this study. The state of sediments in water was iteratively calculated in terms of time, and the action of waves on sedimentary grains was taken into consideration. The simulation only followed fluid dynamics and the patterns of transport, erosion, and deposit of sediments without adding any functional modules that affected the formation and abandonment of river channels or the formation of sand dams and shoals.

Three basic parameters of the delta at different times (i.e., the topset gradient, the number of channel mouths, and delta front rugosity) were monitored during the simulation process. These parameters will be explained and calculated later in this study. When all three parameters of delta morphology became stable, the delta had entered a stable state, and the time when the steady state appeared was at about two thirds of the total simulation time (Figure 5). Therefore, after the simulation of the time duration T, delta morphology was stable and was then analyzed.

Figure 5

Stability analysis of delta numerical simulation. (a) Average topset gradient becomes stable at 0.4 t/T. (b) Average number of channel mouths becomes stable at 0.6 t/T. (c) Delta front rugosity becomes stable at 0.6 t/T. For the sediments of all models, σ=1ϕ, and T is the total time of simulation calculation.

4.1 Comparison of sedimentary simulation of deltas with and without the action of waves

In Simulation Experiment 1, with two different median grain sizes as the inputs, the numerical simulation was carried out with and without wave action. Run IDs A1–A4 in Table 2 correspond to (a)–(d) in Figure 6. The sediment grain sizes shown in Figure 6a and b are both relatively coarse; those shown in Figure 6c and d are also the same, but are comparatively fine. The deltas in Figure 6a and c are river-dominated ones. Those in Figure 6b and d are wave-dominated deltas. It was evident that when there were waves, the number of river channels in the delta was reduced, and in terms of the number of river channels, the wave-dominated delta was more sensitive to changes in sediment grain size than was the river-dominated delta. The actual data of the Mossy Delta in Figure 1a and the Niger Delta in Figure 1b were applied as the simulation parameters in Figure 6c and d. Moreover, the morphologies of the deltas in Figure 6c and d were also close to the morphologies of deltas in Figure 1a and b, which verified the correctness of this wave simulation method.

Figure 6

Comparative simulation of the delta with and without wave action. (a)–(d) Correspond to Run ID A1–A4 in Table 2. The grain size parameters were consistent for (a) and (b) (D₅₀ = 0.5 mm, σ=2ϕ), where (a) Run IDA1 was without wave action, and (b) Run ID A2 was with wave actions. The grain size parameters were consistent in (c) and (d) (D₅₀ = 0.125 mm, σ=2ϕ), where (c) Run ID A3 was without the action of waves, and (d) Run ID A4 was with wave action. The specific parameters of the waves are shown in Table 1.

As shown in Figure 6, once the delta was affected by the waves, not only did the delta morphology change but the distribution characteristics of the river channels changed as well. Because the delta in Figure 6 was the result from the simulation calculation, it could only statically display the distribution characteristics of the river channels. To study the variation of river channels over time during the development phase of the delta, we calculated how the area of active river channels (excluding the abandoned river channels) changed over time (Figure 7). The results showed that the area of the river channels without wave action increased more significantly over time than the area of the river channels with wave action.

Figure 7

Total area of active river channels. The four curves in the figure indicate the change in the area of the river channels over time during the simulation of Run ID A1–A4. Without the action of waves, Run ID A1 and A3 (grey and blue lines) feature a relatively large slope; with the action of waves, Run ID A2 and A4 (red and yellow lines) feature a comparatively small slope

4.2 Comparison of sedimentary simulation of wave-dominated deltas with different sediment grain sizes

In Simulation Experiment 2, based on the nine sediments of different grain sizes input into the simulation system, the delta developed different morphological characteristics under wave action (Figure 8). Figure 8a–i correspond to Run ID B1–B9 in Table 3, from which we directly observed that the shape of the delta and the number of river channels varied significantly with the change in grain size distribution. Moreover, the deltas in Figure 8a, d, and g are quite similar to the three deltas in Figure 2. Therefore, if the causes of changes in delta morphology can be found, the reason for the different morphologies of deltas in Figure 2 can be explained.

Figure 8

Simulation results of sediments with different grain-size distributions under constant wave action. (a)–(i) correspond to Run ID B1–B9 in Table 3. A total of 12,000 steps were involved in the computation, and the simulation time was about 8 years. The specific parameters of the waves are shown in Table 1.

To analyze the effects of sediment grain size on wave-dominated deltas, the morphological features of the deltas were quantified. For the simulation result of each delta, in addition to the three aforementioned basic parameters, i.e., topset gradient, number of channel mouths, and delta front rugosity, the delta shape was also calculated.

To calculate the topset gradient, the coastline in the model area was calculated first. The Opening–Angle method [42] was used in this study to describe the delta coastline. To calculate the coastline, we had to set the threshold value θ of an angle, and here, θ=70°. The coastline can be regarded as the result of smoothing of the 0 m contour line; it was defined as the collection of points whose sum of viewing angles extending to open water, unobstructed by land, exceeded the specified θ. Calculated this way, the coastline can accurately divide the locations of the channel mouths. After determining the coastline, the slope of the line from the apex of the delta to different points on the coastline was measured, and the average was the topset gradient. As the value of ϕ25 of the input sediment increased (R2=0.58), the topset gradient grew smaller (Figure 10a).

Figure 9

Property calculations for the delta simulation results. (a) Delta front rugosity and (b) delta shape; the deltas shown in (a) and (b) are the same, and the sediment lithology could be divided into two types: sandstone 0.1 mm and mudstone <0.064 mm, no waves. The rugosity value of this model ≈ 1.45, A ≈ 0.77.

A channel mouth was defined as a location where more than two channelized cells intersected the shoreline. As the value of ϕ25 of the input sediment increased (R2=0.96), the number of channel mouths decreased (Figure 10b). Corresponding to the data in Table 2, the value of ϕ25 in Figure 8g is 2.917, which is the largest in all the simulations, and there was only one channel mouth.

The value of delta front rugosity was mainly affected by the sand dam at the channel mouth. In this study, the contour line of 1 m underwater was used to describe the delta front rugosity. The contour line was obtained by simply processing the simulated morphological elevation data of the delta. A method similar to Gaussian filtering was used to process the delta front length (red line in Figure 9a), and the filtering window was 40 sampling points long. The delta front rugosity was the ratio of the delta front length to the smoothed delta front length (blue line in Figure 9a). The closer the delta front rugosity got to 1, the smoother the delta front was. If the delta front rugosity was much greater than 1, delta front was rougher. Similarly, as the value of ϕ25 of the input sediment increased (R2=0.71), the delta front rugosity grew smaller (Figure 10c), and the delta front became smoother. Delta shape was quantified with a simple metric:

(9)A=B2L,

where A is the delta shape; B is the delta width, defined as the maximum beach-parallel distance across the delta shoreline; and L is the delta length, defined as the maximum beach-perpendicular distance (Figure 9b). When the A value was greater than 1, the long axis of the delta was parallel to its northern boundary. When the A value equaled 1, the delta was semicircular. Finally, when the A value was smaller than 1, the long axis of the delta was perpendicular to the northern boundary, and the rivers were slender. The simulation experiments showed that as the value of ϕ25 of the input sediment increased (R2=0.46) (Figure 10d), the A value of the delta increased, and the semicircular delta turned into a flat delta with the long axis being parallel to the initial coastline. Run ID B7 had an A value of more than 2, so the delta was approximately arcuate.

Figure 10

Relationship between delta morphological parameters, ϕ25 (a)–(e), and the relationship between the number of channel mouths and Tch¯ (f).

5 Discussion

As detailed in the previous section, we observe a morphological transition in deltas as grain size increases and percent cohesive sediment changes. In this section, we explore how sediment properties create the observed variations in delta morphology when the delta is attacked by the wave.

First of all, the waves affect the delta deposition through the following hydrodynamic processes: (i) waves attacking the river mouth at an angle lead to a deflection of the river jet [43]; they change the jet direction on the river mouth and change the morphology of the mouth bar. (ii) The increase in bottom shear stresses due to waves results in a higher lateral spreading of the jet [44]; an increase in jet spreading produces a sharp decrease of jet velocity, favoring the deposition of the sediments. (iii) The waves increase the maximum shear stress close to the river mouth, maintaining sediment in suspension and transporting it farther away from the river mouth [45]. At last, we must declare that our simulations longshore transport is limited since we consider waves propagating on a small slope bottom. The surf zone is therefore hardly absent, and waves break only near the bar or the beach. So we cannot simulate deltas with strong longshore currents.

5.1 Influence of waves on delta channels

Simulation Experiment 1 consisted of Run ID A1–A4. Figure 6 shows the river-dominated delta (Run ID A1 and Run ID A3) and the wave-dominated delta (Run ID A2 and Run ID A4). Although the river-dominated deltas had different grain-size distributions, they shared a similar river channel distribution (Figure 6a and c), and it was the same case for the two wave-dominated deltas (Figure 6b and d). During the development of the wave-dominated delta, the total area of active river channels grew slowly (Run ID A2 and A4 in Figure 7), whereas the total area of active river channels that were not dominated by waves increased rapidly with time (Run ID A1 and A3 in Figure 7). Figure 7, which shows the changes in the area of active river channels over time, might be too abstract. To more intuitively reflect the influence of wave action on the delta channel, the total duration of river activity was projected on a plane to show the total activity duration of the delta channel (Figure 11). Figure 11a–d are Run IDA1–A4, respectively. The different colors in the figures denote the total activity duration of river channels in the delta during the simulation calculation. The river channels existed in the red area for a long time and existed in the blue area for a short time, but no river passes by the white area. Whether a river channel was active was determined by the water flow in the river channel. If the flow rate was greater than 0.2m/s−1, the river channel was active; if not, the river channel was abandoned. As shown in the figure, in the absence of waves, the river channels frequently bifurcated, and the entire delta was almost filled with river channels (Figure 11a and c). Under wave action, there were fewer but steadier river channels that were only distributed within a small range of the delta (Figure 11b and d). Regardless of the grain size of the sediment, the wave-dominated delta has more stable channels and fewer estuaries than the river-dominated delta.

Figure 11

Time of existence of river channels in the delta. (a)–(d) are time maps concerning the simulation computation of the active rivers Run ID A1–A4 river; (a)–(d) indicate the accumulation of active rivers in terms of time, where red indicates the time of existence of river channels and blue denotes the time period during which the river channels did not exist. The red lines in (a)–(d) represent the delta front.

In our opinion, the direction of the estuaries changes continuously as the delta grows when waves attack the river mouth at low angles, the jet spreading is the dominant process, higher jet spreading favors the deposition of sediments near the mouth, and the formation of a mouth bar. But, the shear stress controlled by the waves and jet has the opposite effect on sediment transport that resuspends and transports it farther away. In this way, the mouth bar has been damaged while being deposited, and the channel cannot branch, making the river channel stable and growing toward the sea. When wave incidence is persistently oblique to the river mouth, low bottom shear stresses favor the deposition of sediments close to the river mouth, likely forming a mouth bar. In particular, waves are perpendicular to the flow at high wave angles, the jet is completely deflected, and bottom shear stresses are low. Littoral drift redirects the sediment downdrift along the adjacent coast. At the river mouth, the side near the beach quickly accepts sediment, forming pronounced subaqueous levees which assume the form of broad shoals. Over time, the natural underwater levees into a delta plain, causing channels that were originally perpendicular to the waves to deflect towards the basin.

The average delta front slope in Run IDA1–A4 was calculated. The delta front slopes in Run IDA1 and A3 without wave action were relatively large, at 0.0942 and 0.0328, respectively. The delta front slopes in Run IDA2 and A4 under wave action were relatively small, at 0.0251 and 0.0061, respectively. This is because the larger the delta slope, the more likely that avulsion would occur. So we think that the slope is one of the reasons why the wave-dominated delta has fewer channels than the river-dominated delta.

5.2 Morphological differences between wave-dominated deltas with different sediment grain sizes

As shown in Simulation Experiment 2, delta morphology changed as the grain-size distribution changed. According to the data statistics, ϕ25 was linearly correlated with the following four delta parameters: topset gradient, number of channel mouths, delta front rugosity, and delta shape (Figure 10a–d).

5.2.1 Relationship between grain size and morphologies of deltas

Changes in delta morphology are associated with sedimentation and erosion. Sedimentation is related to the sedimentation velocity (ws) of the sedimentary grains during migration. The relationship between grain size of the sediment and ws is simple: as the grain size increases, ws increases nonlinearly. Erosion is related to the critical shear stress for erosion (τce), which has a more complicated relationship with the grain size of the sediment. Generally, sediments with a larger grain size (64–500 µm) have a lower value of τce, and sediments with a smaller grain size (≤64 µm) have a higher value of τce [21].

As indicated by the sedimentation simulation results, under wave action, the content of particles with a relatively large grain size in the sediments was the main factor controlling the topset gradient. Sediments with a smaller value of ϕ25 produced a delta with a steeper topset gradient (Figure 10a) because larger particles were less affected by the waves. As a result of the short transport distance caused by a large value of ws, the delta had a relatively large topset gradient. The larger the value of ϕ25, the more significant the action of waves. With a smaller value of ws, more fine-grained sediments were transported to a greater distance, making the topset gradient even smaller.

Why do coarse-grained wave-controlled deltas have more channels than fine-grained wave-controlled deltas? To quantify channel mobility at the top of the delta, we calculated the average time scale Tch¯ of the channels to compare the relationship between different grain-size distributions. Tch¯ denotes the average time period during which the channel remains active in a given position at the top of the delta. Tch¯ was calculated by following the development and subsequent abandonment of river channels. The average time scale Tch¯ in the numerical simulation increased with the increase in ϕ25 (R2 = 0.49) (Figure 10e), which indicated that the delta with a smaller ϕ25 tended to have more frequent avulsion under wave action and more river mouths, whereas the delta with a larger ϕ25 has stable river channels and less river mouths (Figure 10f). This was because the waves had a limited influence on coarse-grained sediments—the delta topset gradient reduction is limited, whereas the steep topset gradient led to an increase in the avulsion frequency of coarse-grained deltas, resulting in an increase in the number of channel mouths in the coarse-grained delta. Furthermore, since the delta formed by finer sediments had a larger τce, that is, the river channels composed of finer sediments had a stronger resistance to erosion, the delta formed by the fine-grained sediments under wave action tended to have more stable river channels.

5.2.2 Different front sedimentation of wave-dominated deltas with different grain sizes

The planar shape of the delta is related to the channel mouth and waves on the coastline. In the simulation experiment, the diversion channel involved the following two situations: (1) avulse to a new location and (2) bifurcate at the mouth bar. Mouth bars are a key factor for delta front rugosity. Delta front rugosity is generated by the difference between the velocities of different parts of the delta. So why do coarse-grained wave-controlled deltas have rugose shape of the delta front with semicircular in planform and fine-grained deltas have smooth delta front with elongate in planform?

In most cases, delta front rugosity corresponds to the sedimentary form of the channel mouth. Although a non-channel mouth area is affected by waves and sedimentation, morphological changes in areas away from the channel mouth are milder than those at the channel mouth. Thus, the morphology in areas lacking channel mouths changes slowly. Near the channel mouth, when the sediment properties forming the delta had a comparatively small ϕ25, a large number of new channels due to avulsion could be found at the delta front.

Figure 12 shows part of Run ID B3, which evolved from a one-river channel (Figure 12a) to an avulsed channel coexisting with the original channel (Figure 12b). In the end, the original channel was abandoned and a new mouth bar originated from the avulsed channel (Figure 12c). In this mode, the delta formed a large number of mouth bars because of a large number of channel mouths, and the mouth bars caused the rugose shape of the delta front.

Figure 12

Formation of avulsed channels dominated by coarse-grained sediments in the delta front (the local part of Run ID B3). (a) step 8160, (b) step 8280 and (c) step 8400.

When the sediment had a relatively large ϕ25, the sedimentary grains became finer, and the delta was more significantly affected by the waves. When the waves spread towards the shore, the mouth bar was formed at the channel mouth (Figure 13a), and there was Channel ①. The underwater distributary channel branched along the mouth bar to develop a new secondary Channel ②, and Channel ① and Channel ② coexisted at this point (Figure 13b). With the passage of time, the mouth bar was transformed by the waves, and the secondary Channel ① near the coastline vanished because of sedimentation (Figure 13c), turning the mouth bar into a natural dyke-dam. Then, as the river jet-deflected towards the direction of the waves (Figure 13d), a new mouth bar was formed and a new secondary Channel ③ was developed (Figure 13e–f). The process repeated until the channel was abandoned (the linearly arranged small black arrows denote secondary Channels ①–③). This process made the delta front even smoother.

Figure 13

Wave-transformed mouth bar. The circled numbers (①–③) in the figure denote the secondary channels of the underwater diversion river, showing the part of Run ID B9. (a) step 3240, (b) step 3480, (c) step 3720, (d) step 3960, (e) step 4200 and (f) step 4440.

In addition, when the sediment had a comparatively small ϕ25, the sedimentary grains were relatively large on the whole, the topset gradient was large, and the delta channel was more likely to have frequent avulsion, making the delta more likely to develop uniformly towards the sea. A large number of channels evenly transported the sediments on the delta coastline and coarse-grained sediments resisted higher bottom shear stress. As a result, the A value of the shape of the delta was relatively small (i.e., close to 1), so the delta resembled a semicircle. When the sediment had a relatively large ϕ25, there was a small number of river channels (Figure 10b), and the mouth bar was frequently destroyed, causing the river channels to stably extend towards the sea. At this time, if more sediments were deposited around the channel mouth, the L value grew larger and the A value was near 1.5; the delta shape approached the shape shown in Figure 8d. If more sediments were carried away by waves, causing most of the sediments to pass around the channel mouth, then the B value of the delta grew larger, the A value of the delta was about 2, and the delta shape was closer to the shape shown in Figure 8g.

6 Conclusions

According to Simulation Experiment 1, a smoother delta front and an unstable mouth bar were formed under wave action. Simulation Experiment 2 showed that because of changes in the grain size, the arcuate delta with a shallow gradient and stable channels turned into a semicircular delta with a steep gradient and many active channels. This morphological change occurred because the nature of the sediments determined the number of channels and the waves reshaped the sediments. Moreover, the ϕ25 of the sediment can serve as an important indicator for studies on deltas under wave action. Based on the results of the numerical simulations and real delta data, the following conclusions are drawn:

1. Affected by waves, the delta front tended to have a small gradient, resulting in the wave-dominated delta to have slender and stable river channels. Less impacted by waves, the river-dominated delta had a relatively large delta front gradient with frequent underwater diversion of channels and avulsed channels.

2. When ϕ25 was small, since the critical shear stress for erosion (τce) and the impact of waves were relatively small, the coarse-grained sediments were input to produce a steep topset gradient. A steeper gradient resulted in a shorter channel time scale (Tch¯) and a larger number of channel mouths. The abundant channel mouths often split and moved, resulting in a nondominant sediment row at each channel and numerous mouth bars, which produced semicircular triangles and the outer boundary of the relatively coarse delta front.

3. When ϕ25 was relatively large, the finer sediments produced a shallow top gradient because of the wave action, resulting in a larger Tch¯ with only a few channel mouths. As a result, the delta had a few channel mouths. Because of the transformation of waves on the mouth bar, natural dyke-dams mainly developed at the channel mouth, destroying the growth of the mouth bar and the division of the river channel. Moreover, coastal sedimentation was caused by wave action. All the factors led to the arcuate delta shape and the smooth delta front.