Horizontal and Vertical Sediment Sorting in Tidal Sand Waves: Modeling the Finite‐Amplitude Stage

The bed of coastal seas displays a large number of rhythmic bed features, of which sand waves are relevant to study from an engineering perspective. Sediments tend to be well sorted over these bed forms, which is so far poorly understood in terms of modeling of finite‐amplitude sand waves. Using Delft3D, we employ bed stratigraphy and consider four different grain‐size classes, which are normally distributed on the phi scale. The standard deviation (sortedness) is then varied, whereas the geometric mean grain size is kept constant. The results show that, typically, the crests of sand waves are coarser than the troughs. Residual flow causes net sedimentation on the leeside of the crest, and, consequently, the general sorting pattern is distorted. Since larger grains experience a larger settling velocity, they are deposited on the upper lee slope, whereas the smaller grains are found on the lower lee slope. Due to sand wave migration, also, the internal structure of the sand wave is revealed, which follows the same sedimentation pattern as the lee slope surface. These results qualitatively agree with sorting patterns observed offshore. The sorting processes lead to longer wavelengths and lower wave heights, as a function of standard deviation. This relates to the dampening effect of suspended sediment transport for fine grains. Finally, it appears that the modeled wave heights fall in the same range as observations in the North Sea. These results are valuable for, for example, predicting the morphological response after engineering activities and determining suitable aggregates for sand extraction.


Introduction
Tidal sand waves are wavy bed patterns found on the bottom of coastal shelf seas, such as the South China Sea (Luan et al., 2010;Zhou et al., 2018), Mediterranean Sea (Santoro et al., 2004), San Francisco Bay (Barnard et al., 2006), and the North Sea . These bed forms have a typical crest-to-crest distance of hundreds of meters and may grow up to several meters in height (Damen et al., 2018a;Van Dijk & Kleinhans, 2005). Sand waves are also mobile, as they display migration rates on the order of 10 m per year in the direction of the dominant tidal current (Besio et al., 2004;Van Santen et al., 2011). They are often superimposed by (mega)ripples, which have heights on the order of centimeters to decimeters (Van Dijk et al., 2008), and are of practical interest since they influence sand transport (Charru et al., 2013). Sediment sorting patterns along sand waves are an often observed phenomenon. Usually, sand wave crests are reported to be coarser and better sorted compared to the troughs (Passchier & Kleinhans, 2005;Stolk, 2000;Svenson et al., 2009;Terwindt, 1971;Van Lancker & Jacobs, 2000), although at some field sites the opposite is observed as well (Anthony & Leth, 2002;. Some monitoring campaigns have also recorded the sediment characteristics on the slopes, and they revealed that the lee slopes are generally composed of finer sediments with respect to the stoss slopes (Cheng et al., 2020). Moreover, field observations have shown that sand wave morphology is strongly related to the characteristics of the sediment, as sand wave heights increase for increasing grain size (Damen et al., 2018a;Ernstsen et al., 2005;Flemming, 2000).
The combination of dimensions and dynamic behavior of sand waves explains they may interfere with various offshore human activities, such as navigation (Dorst et al., 2011) and the construction of pipelines (Németh et al., 2003) and wind farm cables (Roetert et al., 2017). Additionally, the well-sorted crests of sand waves provide a valuable resource for nature-based coastal protection measures (e.g., sand nourishments). Therefore, knowledge on the processes governing sand waves dynamics and dimensions, including its sedimentary structure, is crucial for both the safety of shipping and offshore civil engineering constructions and the formulation of sustainable coastal management strategies.

10.1029/2019JF005430
The initial formation of tidal sand waves has been extensively studied in the past decades by means of linear stability analysis (e.g., Blondeaux & Vittori, 2005;Damveld et al., 2019;Hulscher, 1996). These studies demonstrate how the interaction of an oscillating tidal current with small bed undulations induces a tide-averaged net circulation with near-bed velocities in the direction of the crests of the bed topography and, consequently, sediment transport in the same direction. This mechanism is opposed by gravity, which favors a downslope sediment transport. The competition between these two forces ultimately determines the appearance of the bed forms. Among various extensions of the model by Hulscher (1996), Van der Veen et al. (2006) showed that correcting for grain-size diameters improves the prediction of sand wave occurrence in the North Sea. Later, Van Oyen and Blondeaux (2009a) extended the model with a bimodal sediment mixture (i.e., two grain-size classes) to account for sediment sorting processes. They showed that sediments are indeed redistributed over sand waves, leading to either coarser or finer crests. This mainly depends on two physical mechanisms. First, the balance between hiding/exposure and reduced mobility effects results in coarser (finer) crests for strong (weak) tidal currents. Second, the different sediment fractions exhibit distinct excursion lengths and can thus be transported over a different number of sand waves, counteracting the sorting process. Additionally, Van Oyen and Blondeaux (2009b) showed that the preferred wavelength decreases (increases) due to the presence of a coarser (finer) sediment mixture, which was particularly related to suspended sediment transport. A model-field comparison by Van Oyen et al. (2013) revealed that the sorting processes were fairly well captured by the model. However, a basic limitation of the above-mentioned modeling studies is that they only investigated the initial formation stage, since they are restricted to small-amplitude perturbations. To investigate the nonlinear processes which determine the (equilibrium) dimensions and dynamics, recently several numerical process-based models have been formulated (Campmans et al., 2018;Damveld et al., 2020;Németh et al., 2007;Van den Berg et al., 2012;Van Gerwen et al., 2018). Despite that these models are able to simulate the finite-amplitude stage from a flatbed toward an equilibrium, the modeled sand wave heights are overestimated by several meters. Moreover, nearly all these models neglect sediment sorting processes as they employ a uniform sediment composition. In fact, as of yet there is only one modeling study that investigated the effects of a sediment mixture in the nonlinear sand wave regime , which considered a bimodal mixture similar to Van Oyen and Blondeaux (2009aBlondeaux ( , 2009b.  showed that coarser sediments tend to accumulate at the crests, whereas the troughs stay well mixed. Moreover, contrasting the results from previous studies,  did not find a change in wavelength due to sorting processes. However, this model ignored a few essential processes (e.g., suspended sediment transport and advanced turbulence description) which have shown to be important in modeling the finite-amplitude stage of sand waves (Van Gerwen et al., 2018). Notably, also this study significantly overestimated the equilibrium wave height (up to a factor 2), which could very well be the consequence of the simplified model formulation.
We thus conclude that despite sediment sorting processes over sand waves have been well studied for the initial stage of formation, knowledge on these processes related to finite-amplitude behavior is limited. Furthermore, it is likely that the currently available finite-amplitude models are lacking one or more important processes as they generally overestimate the equilibrium sand wave height. Since it has been shown that suspended sediment transport has a decreasing effect on the equilibrium height (Sterlini et al., 2009;Van Gerwen et al., 2018) and that sediment concentrations in the water column are better approximated by including a mixed sediment composition (Van Oyen & Blondeaux, 2009b), we expect that combining a grain-size mixture, suspended sediment transport, and an advanced turbulence formulation in a numerical model shall lead to a more accurate estimation of sand wave dimensions and dynamics. In addition, representing the bed composition by merely a bimodal mixture seems rather crude, as the bed is much more heterogeneous (e.g., Cheng et al., 2018). Therefore, sand wave models are likely to benefit from a more accurate description of the sediment composition. Altogether, addressing these limitations enables us to expand our understanding of the nonlinear processes governing sediment sorting in relation to finite-amplitude sand waves.
Hence, the goal of this paper is twofold: We aim (i) to understand the effects of a multifractional sediment composition on the nonlinear morphological behavior (wavelength, migration, and equilibrium height) of finite-amplitude sand waves and (ii) to reveal how sediment sorting characteristics (mean grain size and sortedness) differ between the finite-amplitude stage and the (initial) small-amplitude stage of sand wave development. To this end, the modeling approach proposed by Borsje et al. (2013Borsje et al. ( , 2014 and Side view of the water column and stratified bed, illustrating the coordinate system (x, z) with corresponding velocity components (u, w), bed level z b , free surface level ζ, and average depth H 0 . Within each layer (n) the grain-size class (j) is characterized by a grain size d and mass sediment fraction ℱ . The thickness of the layers is indicated by L, where the subscripts a , u , b , tot refer to the active layer, underlayer, base layer, and total layer, respectively. The gray area below denotes the unerodible bedrock layer. (bottom) Deposition and erosion process for a total number of underlayers N = 2. Note that in the deposition example the layers are rearranged due to the new top underlayer ( u 1 → u 2 ) , the lowest one being absorbed by the base layer.
Van Gerwen et al. (2018)-based on the numerical process-based model Delft3D (Lesser et al., 2004)-is used and extended to incorporate multiple sediment fractions and vertical bed stratigraphy. In addition, in this study we investigate symmetrical and asymmetrical forcing conditions, with the latter leading to migrating and asymmetrical sand waves.

Model Formulation
Since we extend earlier work by Borsje et al. (2013Borsje et al. ( , 2014 and Van Gerwen et al. (2018), we restrict ourselves to presenting the details regarding mixed sediment transport and bed evolution. Note that the method described below is readily available in Delft3D (see Deltares, 2012) and that we present the main components which are relevant for this work. For further information about the hydrodynamics we refer the interested reader to one of the above-mentioned references or to the supporting information for a short summary.

Active Layer and Bed Stratification
Let us consider a shallow sea of average depth H 0 , with a sandy seabed consisting of cohesionless sediments on top of an unerodible bedrock layer (Figure 1). Following a 2DV approach, we define a Cartesian coordinate system (x, z) with the x coordinate pointing horizontally and the z coordinate pointing upward, corresponding to the velocity components u and w, respectively. The exchange of sediments between the bed and water column is modeled by means of an extended version of the classical active layer approach (Hirano, 1971;Ribberink, 1987). In the classical approach, it is assumed that the seabed consists of two layers: the active layer (formally the dynamical active layer Church & Haschenburger, 2017) with constant thickness L a and the substrate underneath the active layer. The total available sediment thickness L tot may vary in space and time but is initially spatially uniform. All layers are assumed to be well mixed, and only the sediments in the active layer are instantaneously available for transport.
While the active layer approach is suitable for the initial formation of sand waves, as shown by Van Oyen and Blondeaux (2009aBlondeaux ( , 2009b, the method is unable to describe the vertical sedimentary structure (e.g., the internal history) of a sand wave in the finite-amplitude stage. This is in particular important for migrating sand waves, as the internal structure then gets exposed over time. Therefore, we also apply a layered approach to the passive substrate in order to allow for bed stratification (Deltares, 2012). The substrate is divided into a number N of (well-mixed) underlayers of thickness L (n) u and maximum thickness L u,max and a (well-mixed) base layer. The thickness of the base layer L b follows from the total sediment thickness L tot minus the combined thickness of the active layer and underlayers. Figure 1 illustrates the procedure of deposition and erosion within the layered bed. Deposition leads to a flux of sediments from the active layer toward the underlayers. If the top underlayer reaches its maximum thickness L u,max , a new (top) underlayer is created, and the bottom underlayer merges with the base layer. Conversely, in case of erosion, the active layer will be replenished from the top underlayer. Depending on the amount of erosion, the top underlayer may disappear, and the process will continue with the next underlayer. Note that in case of deposition, this implementation neglects the contribution of the bed load to the sediment flux toward the substrate, which would physically more realistic (Hoey & Ferguson, 1994;Parker, 1991). Nonetheless, this approach has been successfully applied in modeling studies in the riverine (Williams et al., 2016), estuarine (Guerin et al., 2016;Van der Wegen et al., 2011), and coastal (Huisman et al., 2018;Reniers et al., 2013) environment.

Heterogeneous Sediment Transport and Bed Evolution
The seabed is assumed to consist of heterogeneous cohesionless sediments. The bed composition is modeled by introducing a finite number J of sediment fractions with grain size d ( j) and mass fraction ℱ ( ) . Note that these fractions are allowed to vary in time and space, constrained by ∑ J =1 ℱ ( ) = 1. The grain sizes in the sediment mixture can be described by the logarithmic phi scale (Krumbein, 1934): with d ref = 1 mm. The arithmetic and geometric mean grain size (subscript m denotes the mean) are then defined as respectively. The standard deviation d is a measure for the sortedness of the sediment, where d = 0 denotes a perfectly homogeneous bed composition. It is given by Following Van Rijn et al. (2001) and Lesser et al. (2004), bed load sediment transport per fraction S ( ) bed (kg s −1 m −1 ) is calculated according to in which s is a slope correction parameter, calculated by with bs being a user-defined tuning parameter, which is set to 3 following Van Gerwen et al. (2018) and Φ the internal angle of friction of sand (30 • ). Furthermore, is the relative availability of the grain-size fraction at the bed; s the specific density of the sediment, which is equal for all fractions; and u * the effective bed shear velocity. Finally, D ( ) * and T ( j) are the nondimensional bed shear stress and the nondimensional particle diameter, respectively, given by where w denotes the kinematic viscosity of water, w the density of water, g the gravitational acceleration, ′ b the grain-related bed shear stress, and ( ) b,cr the critical bed shear stress based on a parametrization of the classical Shields curve. No bed load transport occurs if b ⩽ ( ) b,cr . Next, the suspended sediment concentration per fraction c ( j) in the water column is calculated by solving an advection-diffusion equation: in which w ( ) s is the friction-dependent sediment settling velocity (see below) and s, x and s, z are the sediment diffusivity coefficients in x and z directions, respectively. These diffusivity coefficients are taken equal to the horizontal and vertical eddy viscosity, respectively, and are thus independent of the grain-size fraction.
Due to the presence of other particles, the settling velocity w ( ) s is potentially reduced. This hindered settling effect (discussed in Van Rijn, 2007) is included as a function of the nonhindered settling velocity w ( ) s,0 and the sediment concentration Here is the reference density, and c tot = ∑ J =1 c ( ) is the sum of the mass sediment concentration c ( j) of all fractions. The nonhindered settling velocity w ( ) s,0 is calculated according to Van Rijn (1993) and reads (valid for d ( j) = 100-1,000 μm) As boundary condition at the water surface (z = ), the vertical diffusive flux is set to zero, that is, The boundary condition at the bed Here, D ( j) and E ( j) are the deposition and erosion rate of the particular sediment fraction, evaluated at the bottom of the lowest computational layer, which is completely above a reference height a above the bed. The reference height a equals 0.01h, with h as the total water depth, and is imposed to mark the transition between bed load and sediment load transport (Van Rijn, 2007). Below, the reference height the concentration is assumed to be equal to the reference concentration, that is, This reference concentration reads

10.1029/2019JF005430
Further details regarding the bed boundary condition and calculation of the concentration profile can be found in the Delft3D user manual (Deltares, 2012). The transport of suspended sediment per fraction S ( ) sus (kg s −1 m −1 ) is calculated by Finally, the sediment continuity equation relates local bed level changes to the divergence of sediment fluxes for each sediment grain class individually and reads in which p is the bed porosity. Furthermore, ℱ ( ) int denotes the mass sediment fraction either in the active layer ℱ ( ) a (during deposition) or in the top underlayer ℱ (1, ) u (during erosion), that is, Summed over all grain-size classes, and using

Spatial Model Setup, Parameter Choices, and Technical Details
The length of the model domain is around 50 km, with a variable grid spacing (Van Gerwen et al., 2018). The central part of the domain comprises 750 evenly distributed grid cells of width Δx = 2 m, which gradually increases to Δx = 1,500 m at the lateral boundaries (see Figure 2a). The finer part of the grid is extended in the direction of the net flow in case of a residual current. In the vertical the model has 60 layers, with a high resolution near the bed (0.05% of the local water depth), gradually decreasing toward the surface. Riemann boundary conditions are imposed at the lateral boundaries to allow outgoing waves to cross the boundary without reflecting back into the domain (Verboom & Slob, 1984). Moreover, the extension of the grid from the area of interest toward the lateral boundaries is sufficiently long to neglect any other boundary-related influences, both for hydrodynamics and morphodynamics. The hydrodynamic time step Δt is 12 s, and a spin-up time of one tidal cycle is performed during which no bed level changes are allowed. To speed up calculations, a morphological acceleration factor (MORFAC; see Ranasinghe et al., 2011, for a discussion) can be used. We have found a MORFAC of 500 to be suitable, since additional runs with a lower MORFAC (250 and 100, presented in the supporting information) did not affect the results in a significant way. As a consequence, one tidal cycle approximates 0.7 years of morphological development. With this MORFAC, simulating 100 years of morphological development takes 5 days on a HPC (high-performance computing) facility. A sensitivity analysis regarding the numerical model setup can be found in Van Gerwen (2016).
Furthermore, the thickness of the active layer L a is set to 0.5 m (see also section 4.2), and a total of 60 underlayers is defined with each a maximum thickness of L u,max = 0.15 m. The thickness of the base layer L b is then defined such that, together with all other layers, the (initial) total thickness of the mobile bed L tot adds up to 20 m.
As explained by Van Gerwen et al. (2018) and Damveld et al. (2020), small variations in model settings may shift the preferred wavelength of the emerging sand wave (fastest growing mode, i.e., FGM). Therefore, we first study the effects of a sediment mixture on the FGM during one tidal cycle, by analyzing the growth and migration rates of a set of small-amplitude sand waves with varying wavelengths, which is further explained in section 3.1. Note that we use MORFAC = 1 in the calculations for the FGM. Subsequently, we use the FGM as a starting point for the finite-amplitude stage, such that we may bypass the initial growth stage-and save valuable computational resources. An example of such an initial bed topography is given in Figure 2a.
The system is forced by the S 2 tidal constituent (angular frequency S2 = 1.45 × 10 −4 s −1 ) in such a way that at the lateral Riemann boundaries a depth-averaged velocity amplitude of U S2 = 0.65 m s −1 is attained. and d = 0.5. The gray columns denote the four grain-size classes used in the numerical experiments. Note that the standard deviation of these four classes is slightly lower due the discrete representation of the PDF. Furthermore, in some simulations a residual current of U S0 = 0.05 m s −1 is superimposed to the tidal forcing. The mean water depth (H 0 = 25m) is equal for all simulations. The roughness of the bed is described by the Chézy coefficient C, which can be calculated by with the roughness height k s as proposed by Whitehouse (2005a, 2005b). For details on the calculation of this ripple predictor, see, for example, Cherlet et al. (2007) and Van Gerwen et al. (2018). Given the settings used in this work, Equation 17 leads to a Chézy coefficient of C = 75m 1∕2 s 1 . Note that we use the mean water depth H 0 to calculate the Chézy coefficient, whereas this should formally be the (space-and time-dependent) total water depth h. However, for the cases considered in this paper, correcting for the total water depth leads to variations in C of only a few percent.
Finally, following Van Oyen et al. (2013), we assume that the sediment composition is normally distributed on the phi-scale. As a consequence, we define a probability density function (PDF) according to For each of the following cases we consider a total of four different grain-size classes ( = 4) with a fixed geometric mean grain size and a varying standard deviation. The mass fractions ℱ ( ) then follow from the PDF, illustrated in Figure 2b for a standard deviation of d = 0.5 and geometric mean grain size of d m = 350 μm. The bin size (Δ ) of each class is set equal to the standard deviation, such that the gray columns denote the discrete PDF of the modeled sediment mixture.   ( j) and mass fractions ℱ ( ) given an increasing standard deviation d . The geometric mean grain size d m for each mixture is 0.35 mm. Note that the additional mixtures used in Figure 8 (based on a different mean grain size) are not included in this overview but that these can easily be calculated following the procedure at the end of section 2.3.

Initial Small-Amplitude Stage
First, we analyze the initial morphodynamic response in terms of the growth and migration rate as a function of the topographic wave number k = 2 ∕ , with wavelength . This is done for a number of sediment mixtures with varying standard deviation. Therefore, the initial wave number k is varied in 10 equal steps from 0.008 m ( = 800m) to 0.063m −1 ( = 100m), after which a fourth-order polynomial is fitted to the resulting growth rates (see below).
Using a complex notation, we write the bed level as follows: where ℜ denotes the real part and A a complex bed amplitude. As a result-and given that small-amplitude perturbations display exponential growth-the growth and migration rate c mig are calculated following Borsje et al. (2014) and read Here, T is the (tidal) period; A 0 and A 1 are the initial and calculated bed amplitude, respectively; and ℑ is the imaginary part. The calculated bed amplitude is determined after one tidal cycle by means of a fast Fourier  , whereas this is 226 to 249 m in the asymmetrical case. This relation between sortedness and preferred wavelength qualitatively agrees to the relation found by Van Oyen and Blondeaux (2009b).
By considering the superposition of an S 2 signal and a residual current ( U S0 = 0.05ms −1 ) , the sand waves display migration. The results show that also the migration rate is affected by the presence of a mixed sediment composition (Figure 3c). Apart from a steady increase in migration with an increasing wave number, the migration rate increases for an increasing standard deviation, as well. The migration rate of the FGM ranges from about 5 m year −1 for a perfectly homogeneous sediment composition to 5.5 m year −1 for a mixture with standard deviation d = 0.5. As one might expect, the symmetrical forcing case displays no migration at all (and thus not shown here).
Next, we analyze the initial sediment sorting processes in the small-amplitude regime. To illustrate this, we present the case with the largest range in grain size ( d = 0.5 ) and focus on the response in terms of the mean grain size and sortedness (i.e., standard deviation). Figure 4 shows the sorting results after one tidal cycle (with MORFAC = 500) for a symmetrical (a and c) and an asymmetrical forcing (b and d). In general, the results show coarser crests and finer troughs, while simultaneously, the crests tend to be better sorted (lower d ) than the troughs, which is qualitatively similar to what was found by . For the other sediment mixtures these observations hold as well, albeit less pronounced. Furthermore, the case with a residual current in the positive x direction displays a small phase shift between the sand wave and mean grain size. Here the coarser sediments tend to accumulate just in front of the crests on the stoss slope of the sand wave. Interestingly, this phase shift is hardly visible for the sortedness, with the highest standard deviation occurring in the troughs, similar to the symmetrical case. This situation arises if the ratios between , where the arrows indicate the direction of the residual current. The top row presents the geometric mean grain size d m , and the bottom row the sortedness (defined as the standard deviation σ d ), where a lighter (darker) shade denotes a higher (lower) degree of sorting. Note that the colormap for panels (a) and (b) denotes a different parameter than the colormap for panels (c) and (d).
the coarse and fine fractions are roughly each other's inverse, which indeed results in a similar sortedness, but different mean grain size.

Development Toward the Equilibrium Stage
Here we first shift our attention to the sorting process of finite-amplitude sand waves, after 20 years ( Figure 5) and 100 years ( Figure 6) of morphological development. Note that, as an addition to the results shown below, a continuous video of the considered run (0 to 100 years) is available in the supporting information.
After 20 years we see that for the symmetrical case the surficial sorting process continues with the same upward coarsening trend as after one tidal cycle. In addition, in the crests region also a vertical sorting is visible, with an increasing upward coarsening. Moreover, the slopes are now better sorted than the crests, whereas the troughs remain well mixed. Figure 5e presents the cumulative sediment fractions in the active layer as defined in Equation 1, where the spatial variation is clearly illustrated. In particular, the finest and coarsest grain-size fractions are affected at this stage. Interestingly, the wavelength of the coarsest fraction is half the wavelength of the finest fraction (and that of the sand wave). This is caused by the fact that in the trough, the largest grains are partly immobile. Hence, the crest-directed transport of coarse grains is lower in the trough than on the slopes, and the coarsest fraction becomes trapped in the trough. Note that this mechanism becomes even more apparent in the equilibrium stage presented further below.
In contrast to the symmetrical case, the sorting process in the asymmetrical case shows some remarkable changes after 20 years, compared to after one tidal cycle. Here we see that the coarsest grain sizes occur on the upper leeside of the crests but that the finest grain sizes occur on the lower lee slope. In other words, the gradient in grain size is much stronger on the lee slope than on the stoss slope. The sortedness displays a strong gradient in the troughs, with a better sorting on the lower lee slope. Moreover, due to the migration (just over 100 m, note the shifted x axis), the internal structure is partly revealed now, and the previously described process consequently creates a vertical sorting as well. The cumulative fractions in the active layer (Figure 5f) confirm these observations and show that particularly the coarsest fraction is much smaller on the lower leeside. Moreover, both the coarse and medium-coarse fraction show a gradual increase in the crest in the positive x direction, with its peak just on the lee side of the crest. For the symmetrical case after 100 years (Figures 6a, 6c, and 6e) the results show the highest grain-size diameter in the sand wave troughs. The vertical sorting structure of the crest reveals a sequential sorting of a relatively coarser, finer, and again a coarser region in upward direction. At the same time, the sortedness still indicates a high degree of sorting on the crests with respect to the troughs.
In the asymmetrical case after 100 years (Figures 6b, 6d, and 6f) roughly three distinct regions with regard to sorting can be observed. Again, also here the troughs display a coarse mean fraction and are well mixed. Conversely, the crests are better sorted but have a relatively coarse mean grain size. In addition, the transport (top) layer on the crests displays a slightly lower mean grain size than the interior. The third region consists of both slopes and the adjacent internal structure and is characterized by a low mean grain size and high degree of sorting. Remarkably, as a result of the horizontal displacement of the sand wave (∼550 m) the internal sorting structure of the sand wave shows hardly any variation in horizontal direction any more, as opposed to the vertical structure.
The large gradients in grain-size fractions on the lee side are the result of the steepness combined with the transition from a mobile to an immobile transport regime in the trough region, which potentially leads to local numerical inaccuracies. Moreover, the absence of lee-face sorting mechanisms may further affect the observed sorting pattern on the lee slope. However, since the observed redistribution process was already visible after 20 years (see the peaks in Figure 5f) and that these gradients gradually increased over time (model output not shown here), we are confident that these numerical aspects do not affect the qualitative process of redistribution. These observations are further illustrated in Figures 7a and 7b, where the sorting characteristics on the crest and in the trough (top layer only) are plotted over time. It is clearly visible that initially, the crests (troughs) become coarser (finer). Once in the troughs the threshold of motion is not exceeded any more for the largest fractions, this process reverses, and a quick increase of both the mean grain size and sortedness occurs (around 30-40 years). Moreover, toward the equilibrium stage, this distribution process stabilizes and results in finer, well-sorted crests and vice versa. It turns out that during this relatively short period of time, the general sorting patterns is formed. Consequently, the sorting time scale is shorter than that of sand wave evolution (see results below).
Next, to show the wave height evolution of different mixtures (d m = 0.35 mm), we present in Figures 7c and 7d the wave height (on a logarithmic scale) as a function of time. In order to isolate the effect of different grain-size mixtures, we fix the initial wavelength of all cases to the FGM of the uniform sediment sample ; see also Figure 3 and Table 1. It appears that for both the symmetrical and asymmetrical cases the mixture only slightly influences the wave height; after 100 years of development the difference between the uniform sample and the one with the highest standard deviation is on the order of decimeters. The final sand wave height is 8.3-8.8 and 7.0-7.4 m for the symmetrical and asymmetrical cases, respectively. It should be noted that, due to the decrease in growth rate, the cases with the highest standard deviation still showed an increase in wave height of a few millimeters per year after 100 years. In addition, Figures 7c and 7d also illustrate the exponential growth during the initial small-amplitude stage, as we assumed for the analysis in section 3.1. In fact, this linear behavior continues until a sand wave height of around 4.5 m is reached, after which nonlinear effects become increasingly dominant and the exponential growth is eventually damped.  Figures 4-6). Note that the actual standard deviation is slightly lower, as illustrated in Figure 2b. In panels (c) and (d) the wave height development over time for a range of different standard deviations with d m = 0.35 mm. The vertical axis is plotted on a logarithmic scale to emphasize the linear behavior during the early growth stage.
From the initial formation process it is known that in the currently employed shallow water model, a uniform grain size of around 200-250 μm leads to negative growth rates, and thus inherently a flatbed (Borsje et al., 2014). This motivates us to study the sensitivity of the equilibrium wave height to the geometric mean grain size. Again, we fix the initial wavelength to that of the FGM for a uniform mixture, and we consider three different standard deviations, namely, d = 0, 0.5, and 0.75. Figure 8a shows a substantial decrease in sand wave height with decreasing mean grain size. In the case with a standard deviation of d = 0.5, the range in wave height after 125 years is 4.6-7.4 m for a geometric mean grain size of 0.25-0.40 mm, respectively. In particular, for the smallest grain sizes the decrease in wave height is strongest. As expected, due to the increasingly larger dampening effect of the (fine) sediment fractions in suspension (Borsje et al., 2014) (to be discussed further below), the sand wave for the mean grain-size case of 0.20 mm decays toward a flatbed. The deformations which are visible between 75 and 100 years of development for the coarsest grain-size cases are the result of mobility effects, as at this approximate (trough) depth, the critical shear stress is not exceeded any more for the largest fractions. Hence, the troughs are not eroding any more, and the overall growth rate initially decreases, after which the crest growth increases slightly to compensate for this effect. Note that the different depths at which these deformations occur is due the differences in critical shear stress between the fractions. Finally, similar to the results from Figure  7, an increasing standard deviation leads to a further decrease of the sand wave height (Figure 8b), which is further discussed in section 4.1.1.

Discussion
We extended the nonlinear sand wave model of Van Gerwen et al. (2018) by including multiple heterogeneous sediment mixtures and bed stratification, which allowed us to investigate sediment sorting processes over finite-amplitude sand waves. We showed that during the earlier growth stages, the troughs of sand waves tend to be finer than the crests, whereas the opposite holds in the equilibrium stage. Next, increasing the degree of sorting leads to an increasing reduction of the modeled sand wave height. Below, we discuss field observations, the physical interpretation of the results, and the limitations of this paper. Damen et al. (2018a) present a large-scale analysis of aggregated data sets on the relation between various environmental parameters and characteristics of sand waves. The data consist of almost 10,000 samples from the Dutch Continental Shelf and has a resolution of one data point per square kilometer. For a model-field comparison we use data on the observed wave height and sediment diameter, which are plotted in Figure 8b. A least squares regression line is plotted through the data to indicate the relation between the two variables. As already pointed out by Damen et al. (2018a), the (positive) slope between sand wave height and grain size is steeper for smaller grain sizes. To emphasize this, we divided the data in two bins, namely, d m < 0.35 mm and d m ⩾ 0.35 mm.

Sand Wave Height
In Figure 8b also, the modeled wave heights after 125 years (from Figure 8a) are shown. Although close to the upper limit of the field data, the trend of the model results agrees well with the field observations. Particularly, the increasing wave height with respect to an increasing grain size is well captured by the model. Moreover, the other large black markers in Figure 8b clearly indicate a decrease in wave height for an increasing sortedness. This effect is strongest for the most fine and coarse mean grain sizes. Similar to the deformations in Figure 8a, the seemingly deviant result from the case with d m = 0.4 mm and d = 0.75 is due to the transition from a mobile to an immobile bed.
One should note that each data point represents a different set of field conditions, whereas the model setup only distinguishes in sediment characteristics. Previous research has revealed that tidal asymmetry and waves also lead to significantly lower sand waves (Campmans et al., 2018;Sterlini et al., 2009;Van Gerwen et al., 2018). Moreover, other environmental parameters, such as water depth, current velocity, and the Chézy coefficient (which is in turn a function of depth and sediment diameter), are likely to affect the wave height, too, while these parameters were kept constant here.
In general, the model results showed a reduction in wave height as a result of the presence of a sediment mixture. In particular the smaller fractions in the mixture are responsible for this reduction, as smaller sediment grain sizes generally lead to lower wave heights (e.g., Figure 8a). This is caused by the dampening mechanism of suspended sediment transport, which is the result of the highest suspended sediment concentrations occurring downstream of the sand wave crest (Borsje et al., 2014). As this phase difference is present both during flood and ebb, the net sediment flux is away from the crests. Moreover, Borsje et al. (2014) showed that this effect is inherently larger for smaller grains, ultimately leading to sand wave decay for small-sized sediment mixtures. The distribution of the different grain-size fractions in the active layer (Figures 5e, 5f, 6e, and 6f) corroborate this, since they clearly show the accumulation of the finest fraction in the trough region of the sand waves, highlighting the dampening mechanism of fine grains.
As already shown by Van Gerwen et al. (2018), tidal asymmetry leads to a further reduction in wave height. Due to the presence of a residual current, the tide-averaged circulation is distorted, and hence, convergence of sediments does not occur exactly above the crests. Here, this observation is confirmed by the modeled sorting pattern. In the symmetrical case, the coarse sediments clearly accumulate on the crests of the sand wave (Figure 5a), whereas this accumulation is on the leeside slope of the crests in the asymmetrical case (Figure 5b).

Sediment Sorting
Here it is not our aim to perform a site-by-site comparison of field data and the modeled sorting results. From the field sites where extensive grain-size measurements are available, other environmental data are generally not detailed enough to provide a well estimate for, for example, the local hydrodynamic conditions. Hence, calibrating the model to a specific site may lead to large uncertainties in the magnitude of the modeled grain sizes. Instead, we look for field evidence that is able to confirm the general sorting patterns presented in this work.  and Van Oyen et al. (2013) present a comprehensive overview of surface sorting patterns over sand waves in the North Sea. Seven out of the 10 reported field sites in these studies are characterized by a coarser crest, compared to the troughs. Although our model results show a coarsening of the trough in the equilibrium stage (due to the influence of the critical shear stress; see further below), the sorting process during the earlier growth stages agrees well with these observations. Both in the symmetrical and asymmetrical cases, the overall trend in the model is a coarser crest with respect to the troughs. Moreover, the rate of sorting in the field indicates a high sortedness (low standard deviation) on the crests, which again agrees well with the model results. However, the eventual mismatch between the model and the field observations suggests that some processes are missing or not adequately included.
Of the three reported field sites with finer sediments at the crests, two were located superimposed on larger sand banks , whereas the other was influenced by a strong, episodically ocean current (Anthony & Leth, 2002). These background effects hinder the assessment of the responsible mechanism for this observed crest fining, as they are likely to influence sediment transport patterns in these areas to a significant extent. Despite, our model also showed that finer crests with respect to the troughs could occur, which is the result of the coarsest fraction becoming immobile. It turns out that, when the trough depth increases, the critical shear stress in the troughs is no longer exceeded for increasingly longer periods of the tidal cycle. This is clearly visible in Figure 6e, where the coarsest fraction makes up more than 60% of the sediments in the trough, whereas the medium-coarse fraction (which is still mobile) is almost fully eroded. Note that the finest fractions are still present in the troughs due to the dampening contribution of suspended sediment transport. Nonetheless, based on the available data from these field surveys, it is unclear whether the observed trough coarsening was the result of mobility effects or other processes. In this respect, we recommend to further study the role of the threshold of motion in this matter by comparing alternative transport formulations which do not consider this particular process (e.g., Wilcock & Crowe, 2003).
In the model by Van Oyen and Blondeaux (2009b) and Van Oyen et al. (2013) finer crests were the result of weak tidal currents which favor the transport of fine materials. Indeed, this is similar to what our results suggest, as either weaker currents or coarser sediments lead to lower transport rates for these coarsest fractions.
In fact, our model also revealed (not shown here) that the initial redistribution of bed materials-instead of the showed equilibrium situation-may lead to finer crests if initially the coarsest fraction is (partly) immobile. Strictly speaking, this is even better comparable to the results of Van Oyen and Blondeaux (2009b) and Van Oyen et al. (2013) since they only studied the initial formation process.
However, it is questionable if this particular case of crest fining is the explaining mechanism for the mentioned field observations, because of the large uncertainties involved in the data as discussed above. Moreover, the coarsening of the trough found in the model is more likely to be analogous to the presence of unerodible (armored) layers in the bed, as for instance observed in the English Channel (Le Bot & Trentesaux, 2004) and modeled by Porcile et al. (2017).

10.1029/2019JF005430
Other than crest/trough sorting, observations on sorting patterns over the slopes are more scarce. Cheng et al. (2020) present a detailed field campaign where they found that the stoss slopes of sand waves are coarser than the leeside slopes. In case of a superimposed residual current our model shows comparable results, but only for the initial response (Figure 4b). During the later stages, the peak of mean grain size occurs on the upper lee slopes, whereas the lower lee slope has a finer mean grain size (Figure 5b). Since the study by Cheng et al. (2020) did not report on hydrodynamic conditions, it is unknown what caused this contradiction. However, an explanation for the model result might be found in the fluvial environment. The governing mechanisms for dune formation in rivers and tidal sand waves show quite some similarities (Hulscher & Dohmen-Janssen, 2005). In case of a superimposed residual current, both sorting mechanisms show an insightful resemblance (e.g., Kleinhans, 2004, and references therein). Grains with the largest diameter experience the highest settling velocity (see Equation 9). Therefore, these larger grains are the first to be deposited when the flow velocity decreases, which is on the upper lee slope. Conversely, the finest grains are deposited on the lower lee slopes. Note that if the lee slope angle is high enough, possible other mechanisms which are excluded in the model (e.g., avalanching) may also affect the sorting process. Moreover, it is recommended to further study the role of the slope correction parameter s (Equation 4) on the sorting processes in the slope region, similar to what was done by Wang et al. (2019) for wavelength and migration in the initial formation stage.
Finally, as a result of migration the model results also reveal the internal structure of the sand waves. This structure follows the sedimentation pattern described above, where the coarsest grains accumulate in the upper part of the sand waves and the finer grains in the base of the sand wave. In the equilibrium stage (Figure 6b), the active layer at the stoss slope displays a finer mean grain-size pattern than after 20 years of development ( Figure 5b). Due to the horizontal displacement of the sand wave, the finer grains originating from earlier deposits become mixed in the top layer, leading to a fining of the stoss slope. In contrast to the lack of field evidence for sand waves, a study by Koop et al. (2019) into migrating mega ripples did reveal a similar sorting pattern where both the leeside and stoss side are finer than the crest and trough. However, all these field data were obtained by box corer and multicorer, which only describe the upper part (∼20 cm) of the bed, so that the internal structure cannot be assessed. Hence, there is a need for field evidence of deeper sediment layers (using, e.g., vibrocorers) that uncovers the internal structure of sand waves in terms of mean grain size. Alternatively, ancient dune deposits in terrestrial environments (e.g., the Belgian Brussels Sands Houthuys, 2011) can also provide insight regarding the internal sorting of sand waves.

Assumptions and Limitations
The majority of modeling studies into grain-size sorting over marine bed forms considered the hiding-exposure mechanism (Foti & Blondeaux, 1995;Van Oyen & Blondeaux, 2009a;Van Oyen et al., 2011;Roos, Wemmenhove, et al., 2007;Walgreen et al., 2004). It assumes that the transport of finer grains is hindered, as they are protected by movement from the surrounding larger grains. Conversely, larger grains are more exposed and therefore transported more easily. Our model did not include this mechanism, so that theoretically the transport rate of fine grains is overestimated and that of the coarse grains underestimated.  pointed out that excluding the hiding-exposure formulation does lead to comparable results for the initial stage of sand wave formation and that only the rate of the initial sorting process was enhanced. Moreover, we showed that our model, by only taking into account sediment mobility effects, is able to represent the observed sorting patterns in sand wave areas, so that the hiding-exposure formulation is thus not essential here. Therefore, it is likely that including such a correction to the sediment transport will lead to differences in both the sorting rate and quantitative grain sizes but that qualitatively the observed sorting patterns will not change. Nevertheless, future efforts focusing on sorting processes in the finite-amplitude stage should include further study on the quantitative effects of the hiding-exposure mechanism.
Under certain conditions, the active layer model may become ill posed (Chavarrías et al., 2018;Ribberink, 1987;Stecca et al., 2014). In these ill-posed situations, short-wavelength perturbations may be triggered (by for instance numerical noise) and grow until a point where the mathematical problem becomes well posed again. In a practical sense, these oscillations lead to changes of the stratified bed, which are physically not valid. According to Chavarrías et al. (2018), it is typical for ill-posed problems that finer grids result in larger errors since the erroneous perturbations have more space to grow, whereas solutions of well-posed problems converge with decreasing grid size. For the current model, it turns out that both the bed level and the sorting properties are not affected by a decreasing grid size (and time step), such that we are confident that the presented results are physically correct. As an alternative to our employed method, recent advancements have resulted in approaches (e.g., the SILKE model Chavarrías et al., 2019) that ensure the well posedness of the active layer concept.
The thickness of the active layer is assumed to be in the range between a few grain-size diameters  and the height of the superimposed bed forms (Van Oyen & Blondeaux, 2009b). Here we have set this thickness to a spatially uniform value of 0.5 m, which is generally larger than the bed forms over sand waves, although megaripples on sand wave crests can have wave heights of several tens of centimeters. Still, since in sand wave troughs ripples are usually much smaller or even absent , it is likely that the modeled active layer thickness is overestimated in the model. Importantly, this thickness is related to the time scale of the sorting process , where a larger thickness increases the time scale. Here, the results (see also the videos in the supporting information) show that the sorting time scale is already much shorter than the time scale of sand wave evolution. An even shorter sorting time scale may lead to a faster coarsening of the trough and, consequently, hinder sand wave development. Therefore, we performed additional model runs (presented in the supporting information) with an active layer thickness of 0.25 and 0.1 m. It turns out that the sorting pattern is unaffected by this decrease of the active layer thickness.
To determine the effect of sorting processes on the wave height, we started each simulation with a fixed initial wavelength. However, our results also showed that the preferred wavelength (FGM) is affected by changing sediment characteristics (see also Van Oyen et al., 2013, andWang et al., 2019). Particularly, the wavelength increases for increasing standard deviation (lower sortedness) and for decreasing mean grain size. Since the equilibrium height is potentially affected by a changing wavelength (Blondeaux & Vittori, 2016;Damveld et al., 2020;Van Gerwen et al., 2018), one should realize that the presented results for the wave heights are of a particular mode, rather than that of the FGM of that case. Moreover, section 4.1.1 already discussed many other variables that affect the wave height. This complexity strengthens our choice for a fixed mode, which conveniently allows to isolate the relation between sediment sorting and finite-amplitude behavior.
We assumed uniform conditions for the initial sediment composition, that is, a well-mixed cohesionless sandy bed. However, in reality the bed consists of a wide range of sediments including mud fractions  and displays a large spatial variation, such as the above-mentioned bed armoring (Le Bot & Trentesaux, 2004) and incidental clay layers (Allen, 1982). Besides, these spatial variations are also expressed through the nonuniform distribution of benthic organisms over sand waves (Damveld et al., , 2020, and these organisms are able to influence sediment transport processes to a large extent (Malarkey et al., 2015;Widdows & Brinsley, 2002). Furthermore, also the Chézy roughness coefficient (Equation 17) was assumed spatially uniform in this study, whereas it is known that seabed roughness varies over sand waves . Moreover, as follows from the results, skin roughness is also nonuniform over sand waves. In turn, these spatial variabilities in seabed roughness may significantly affect sediment transport processes. Hence, including these nonuniformities in the present numerical model may even further explain local variations in the observed sorting patterns.
Finally, the Delft3D model (Lesser et al., 2004) is based on the shallow water assumption; that is, the horizontal length scale (wavelength of sand waves) is much larger than the vertical length scale (water depth), such that vertical momentum equation can be reduced to an expression of hydrostatic pressure. Thus, vertical accelerations are assumed to be small compared to the gravitational acceleration and can therefore be neglected. However, it is questionable whether these vertical accelerations remain small for equilibrium sand waves, and in particular near the steep (lee) slope. Hence, it is recommended to study the differences between nonhydrostatic and hydrostatic flow over equilibrium sand waves, by, for instance, adopting a similar approach as used by Lefebvre et al. (2014).

Practical Implications
The current model can be applied for predicting site-specific characteristics of sand waves and, ultimately, determining the morphological effect of human-induced disturbances in the marine (sand wave) environment. Additionally, including other known processes that influence nonlinear sand wave morphology (e.g., wind effects (Idier et al., 2002)) would even further increase the predictive capacity. Next, this model can support benthic habitat mapping, as sediment characteristics are a good indicator for the presence of 10.1029/2019JF005430 benthic organisms (Reiss et al., 2014), thereby providing better knowledge for sustainable ecosystem management. Finally, this model can be used for determining suitable aggregates for extraction, information which is valuable in a marine spatial planning context.

Conclusions
We presented a morphodynamic model which is able to describe sediment sorting patterns in finiteamplitude sand waves in a more sophisticated way than previously applied in literature. In particular, we used a multifractional approach to describe the sediment composition and applied bed stratigraphy to uncover the internal sorting structure of the sand waves.
The model results showed that, in general, the troughs of the sand wave hold finer sediments than the crests, which corresponds well to field observations on sediment sorting. This general pattern is distorted due to tidal asymmetry, which leads to an accumulation of coarser sediments on the upper lee slope, and finer sediments on the lower lee slope. Due to sand wave migration, this sedimentation pattern consequently also determines the internal structure of the sand wave. Furthermore, since the coarsest fractions are the ones to first experience immobility by not exceeding the critical shear stress, the sand wave troughs eventually become coarser toward the equilibrium stage, and an unerodible layer is formed.
The morphological response to an increasingly less-sorted sediment mixture results in longer wavelengths and lower wave heights. Moreover, by varying the geometric mean grain size of several heterogeneous sediment mixtures, the model provides estimates of wave heights that are in the range of observed sand wave heights in the North Sea.

Data Availability Statement
Data are available through Damen et al. (2018b).