Explaining the formation of sedimentary structures under antidunes using a 2D width-averaged numerical model

: Explaining the formation of sedimentary structures under

Observations of the hydraulic conditions during the formation of geological records are usually not available, as the depositional processes often are in the past. Eye-witnesses to large floods may exist (Duller et al., 2008), but their observations are limited to the shape of the water surface. The sediment processes taking place at the bed are most often hidden under sediment-laden water in a field situation (Alexander & Fielding, 1997). Near-bed processes can be observed in laboratory channels with glass walls (Alexander et al., 2001). There is, however, a limit to how detailed flume studies can describe complex sediment transport processes. Even in a controlled laboratory environment, it is difficult to measure velocity, pressure, turbulence and shear stress near a moving bed or the water surface.
However, all these parameters can be computed numerically in a technology called Computational Fluid Dynamics (CFD) (Vellinga et al., 2018). A CFD model that computes the sedimentary structures just like a flume study will give much more detailed information about the formation process. The model can visualize temporal and spatial variations in grain size distributions of the bed. If the colour pattern identifying sediment structures in the field can be correlated with the computed grain size distribution (Bridge & Best, 1997;Hofstra et al., 2018), then the computer model can predict the formation of the sedimentary structures. This enables a more comprehensive and detailed analysis of how the sedimentary structures are formed. This is the main motivation for the present study and its novelty compared with previous research.
CFD methods have been used over a number of years in hydraulic engineering (van Rijn, 1986, Celik & Rodi, 1988Olsen & Melaaen, 1993;Savage & Johnson, 2001;Zinke et al., 2011). The focus of these studies has often been to test how well the numerical model can replicate measured velocity profiles or bed elevation changes. The purpose was to assist engineering design of hydraulic structures.
Numerical models have also been used to explain physical phenomena related to geomorphology. Booker et al. (2001) explained shear stress distributions in a natural river with pool-riffle sequences by using a CFD model. Olsen (2020) used a CFD model to explain the avulsion process in a braided river.
Detailed figures were presented, showing the effect of water elevation changes and secondary currents over time. This gave a detailed description of the physical processes occurring in the avulsion, with more parameters than obtained in earlier physical model studies. The present study has a similar goal: To use the CFD technology to describe the formation of the sedimentary structures in more detail.
Numerical methods have earlier been used to model aquatic sand dunes in rivers and channels with subcritical flow. Naqshband et al. (2015) and Warmink et al. (2014) computed formation of dunes by a 1D numerical model solving the shallow-water equations. Dore et al. (2016), and Giri & Shimizu (2006) used a more accurate 2D width-averaged CFD model based on the Navier-Stokes equations. Fully 3D numerical modelling of aquatic dunes in rivers and laboratory flumes has also been done by several researchers (Rüther et al.,2008;Goll, 2017;Nabi et al., 2013;Mewis, 2016;Kopmann, 2020) However, none of these studies computed the detailed sedimentary structures of the deposits.
Antidunes are more difficult to model than dunes, both in the laboratory and with numerical models. This is due to the more complex geometry of the water level. Early laboratory studies of antidunes include Kennedy (1960Kennedy ( , 1963. He also carried out mathematical/theoretical analysis of these bed forms, giving formulas for characteristic parameters, for example the wavelength. Later flume studies with antidunes include Alexander et al. (2001), Recking et al. (2009) andCartigny et al. (2014). Núñez-González & Martín-Vide (2011) showed that it was possible to obtain stable downstream migrating trains of antidunes with regular size persisting over a longer time period.
Numerical modelling of antidunes is less common. Computation of water and sediments is not straightforward, due to complex interaction between the two phases during erosion and deposition.
While a horizontal water surface may be assumed when modelling aquatic dunes, antidunes requires an advanced numerical model to find the vertical location of the free surface. This surface will have a more irregular shape, especially if wave breaking is to be included. Vellinga et al. (2018) computed cyclic steps using the commercial CFD program Flow3D. The air-water interface was modelled with a twophase volume of fluid approach, and a fixed 2D width-averaged grid. The cells could then be filled with the three phases: air, water or sediments. A three-phase numerical model was then required, which is reasonably complex. Olsen (2016Olsen ( , 2017) used a different approach where an adaptive grid tracked the movement in both the free surface and the bed. This simplified the numerical algorithms for solving the Navier-Stokes equations, but increased the complexity of the boundary conditions of the bed and the surface, and the grid regeneration. The latter approach is used in the present study.
The two earlier studies of downstream-migrating antidunes (Olsen, 2016(Olsen, , 2017 replicated hydraulic conditions from the laboratory studies of Kennedy (1960) and Núñez-González (2012). The flumes in these two projects had different bed slope and bedform height. In the present study, hydraulic input parameters similar to the laboratory study of Núñez-González (2012) was chosen. Here, the bedforms were higher than in the experiments by Kennedy (1960), giving a more pronounced effect on the stratification of the sediments.
Earlier numerical studies computing sedimentary structures from dunes, antidunes and cyclic steps only used one sediment fraction (Vellinga et al., 2018). The novelty of the present article lies in the use of a new method to compute the vertical stratification of sediment deposits. The bed layer is divided into 1000 vertical layers, each containing information about sediment fractions of five different sediment sizes. The vertical size of each layer changes over time as a function of sediment deposition and erosion. Thereby, the stratification is recorded. Graphics show in detail how the sediment structure is formed over time, including location and magnitude of beds, laminae and bounding surfaces, together with water depth, velocity, turbulence and bed shear stress. This is used for a detailed explanation of the formation of sedimentary structures observed in the field.
In the following, the paper first explains the numerical model for water flow, sediment transport and formation of the sediment stratification. Then information about the numerical setup for the antidune case is given, and the results are presented with figures and explanation of the sedimentation processes. Finally, a discussion and conclusions are presented.
Description of the numerical model The numerical model contains multiple grids and algorithms. This is explained in two sections: First, the algorithms for computing the water flow and sediment transport are given. Then the new algorithms for the sediment stratification at the bed are described.

Numerical model for water and sediment flow
A two-dimensional width-averaged grid was used to compute the water and sediment flow over the antidunes. The grid was non-orthogonal and adaptive, following both the free surface and the bed level with a constant number of cells over the depth. Fig. 1 shows a longitudinal profile of a typical section of the grid modelling antidunes.
The numerical model computed the water velocities from the Reynolds-averaged Navier-Stokes equations. The k-ϵ model (Launder & Spalding, 1974) was used to compute the turbulence.
The pressure was computed with the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) (Patankar, 1980). A finite volume method was used for discretization of all equations. A bed roughness value, k s , was required for the wall laws (Schlichting, 1979) at the bed. This was used as a boundary condition in solving the Navier-Stokes equations: (1) The velocity, U(y), is a function of the vertical distance from the bed, y, and the shear velocity, u * .
The is a constant equal to 0.4. All parameters are given in the SI system, so k s and y are in meters.
The following formula was used to determine the bed roughness (van Rijn, 1982): The parameter d 90 is the sediment size where 90% is finer, based on the grain size distribution of the top sediment bedlayer. The bed shear stress was computed from the turbulent kinetic energy, k, close to the bed (Rodi, 1980): The k value was found from the k-ϵ turbulence model, where the parameter c is a constant equal to 0.09.
The free surface location was computed with the CGA (Continuity Gravity Adaptive grid) method (Olsen, 2015). The algorithm calculated the water surface elevation changes from the continuity equation in the cells closest to the surface. Normally, the SIMPLE method uses the water continuity defect to compute the pressure. In the CGA method, the water continuity defect in the top cells is instead used to move the water surface vertically. The adaptive grid followed the water surface, as shown in Fig. 1. The sediment transport was computed from multiple convection-diffusion equations, one for each sediment size. The Engelund-Hansen (Engelund & Hansen, 1967) sediment transport formula was used as boundary condition for the bed: (4) The sediment transport q s [kg/s/m] for sediments with diameter d is given as a function of the bed shear stress, , depth-averaged velocity, U and gravity, g. The density, s , of the sediment particles was set to 2650 kg/m 3 , and the water density, , was 1000 kg/m 3 .
The density of the water/sediment mixture was assumed to be constant, meaning no effect of temperature was included. Also, the sediment concentration was assumed to be so small that it did not affect the fluid density. The water content of the bed was set to 50%. The bed movements were computed from the change in the sediment transport capacity in the flow direction.
Other boundary condition used in the numerical model was a constant value of the inflowing water discharge and sediment concentration. Zero gradient boundary conditions were used at the outflow section. A zero gradient boundary condition was also used for the horizontal velocity and turbulence at the water surface.
The solution method and the numerical algorithms described here had previously been tested on several different flow problems (Wilson et al., 2003). The main reason for using this model in the present study was that its ability to compute the formation of antidunes. Olsen (2016Olsen ( , 2017 showed that the computed height, length and the celerity of the antidunes compared well with laboratory experiments. These earlier studies looked at an equilibrium situation, without net sediment accumulation or erosion at the bed. Additionally, only one sediment size was modelled. Therefore, the sedimentary structure based on a variation in the grain size distribution of the bed could not be computed.

Numerical model for sediment stratification
Observations of sedimentary structures in outcrop studies are based on interpretation from eye observations or photographs. An example is given in Fig. 2. Bedding and laminae are then seen as darker or brighter lines. In the present study, it is assumed that similar grain sizes will have the same colour or greyscale (Bridge & Best, 1997;Hofstra et al., 2018). Lines of particles with similar particle diameter will then form the sedimentary structures. The main novelty of numerical model presented here is its ability to compute the grain size distribution of the bed with a high vertical resolution. This means that the numerical model can directly compute the sediment stratification over the depth.
The sediment stratification algorithms were based on a very fine mesh in the vertical direction.
The present study used 1000 vertical bed layers. Since the grid had 700 cells in the streamwise direction, the number of bed layer cells resolving the sedimentary structures was 700 000. The layers were initially distributed uniformly in the vertical direction. Each layer contained the information about the sediment fraction of each sediment size. The thickness of each of the 700 000 layers could also be varied individually.
The layers changed as erosion and deposition occurred. If erosion was computed at one time step, the top bed layer would be reduced in size or completely removed. If instead deposition occurred, a new layer would be formed. The number of vertical layers was finite (1000) and one model run included around 1 million time steps. The formation of a new top layer meant that two of the layers below had to be merged. The question was which layers? An algorithm was used to search all layers to find the lowest difference in the grain size distribution. This was computed from the following formula: The index k denotes the vertical layer number, varying from 1 to the 999 layer interfaces. The index j denotes the sediment size, from 1 to 5 sediment sizes. The sediment fraction, f, is a floating point number that varies between 0.0 and 1.0. In Equation 5, f has subscripts with sediment size j, and layer number k.
The two layers with the lowest G value in the vertical column were then joined. This new layer had a vertical size equal to the sum of the two original layers. Merging the two layers gave room for the one new depositing layer in the stack of 1000 layers. Since the joining operation was done in the layers with the lowest fractional differences, the stack of layers formed large cell heights in areas of more uniform grain size distributions and smaller cells where the gradients in f were steeper. An illustration is given in Fig. 3, showing a small section of the bed sediments.
The numerical model also included an algorithm for sand slides. It is possible to compute small avalanches with complex geotechnical methods (Olsen & Haun, 2020), but in the present study a simpler option was chosen. The bed elevation of two neighbouring cells were compared to check if the bed slope exceeded the angle of repose for the particles. This was set to 31 degrees in the present study. If the local bed angle exceeded the angle of repose, the surplus sediment from the higher cell was moved to the lower cell. This was done by eroding a part of the top layer from the upper cell and depositing it in the lower cell as a new layer. The same algorithms was used as in the computation of erosion and deposition by moving water. This meant sediment continuity for the different fractions was preserved. Also, stratification due to sediment slides was modelled, as shown later. The resulting bed slope became equal to the angle of repose. This procedure was repeated for all cells after each change in the bed elevation at the end of a time step.

Modelling stationary and downstream migrating antidunes
The present study used the input data and geometry from the S2 flume study of Núñez-González  (Olsen, 2016). The grid resolving the water flow field had a cell size of 1 cm in the streamwise direction and initially 0.4 cm in the vertical direction. This meant 25 cells over the water depth.
A very short time step of 0.0005 seconds was needed for stability reasons. The computational (wall clock) time on a 3.4 GHz PC was typically four weeks for one simulation with 1000 bed layers, modelling deposition over 5000 seconds, or 1 million time steps. The choice of 5000 seconds was based on the need to generate a significant thickness of the bed deposits but also not to have too excessive computational times.
The original bed sediments from the study of Núñez-González (2012) had a uniform diameter of 1.48 mm. The goal of the current research was to model sedimentary structures, requiring a non-uniform sediment distribution with multiple sizes. The five sediment sizes had diameters 0.5, 1, 1.5, 2 and 2.5 mm. An equal inflow of each of the sediment fractions was used, so that the average sediment diameter would be the same as in the laboratory study. The physical model study used a sediment inflow concentration of 5.8 g/l. This number was found to give an equilibrium bed elevation, both in the physical model and in the previous numerical modelling study (Olsen, 2016). The interest in the present study was to investigate sediment deposition. Therefore, the inflowing sediment concentration was increased to a higher value of 9.3 g/l. This caused a reasonable amount of bed aggradation, while the hydraulic parameters did not deviate substantially from the equilibrium values.
The average water depth and the water velocity remained reasonably constant throughout the simulation. The measured average water depth from the physical model study of Núñez-González (2012) was 0.109 meters. This was 1.7 cm higher than the results given by the numerical model, which was 0.092 meters. The depth values from the numerical model had been averaged over the last 4000 seconds. There were fluctuations in the computed water depths, and a value over 0.109 meters was often reached. The length of the antidunes in the physical model study was 0.52 meters, 8% longer than the average value predicted by the numerical model, 0.48 meters. The observed bed form height from the laboratory was 0.065 meters, while the average value from the numerical model was lower, 0.056 meters.
The increase in bed level at the upstream and downstream parts of the flume was also computed over time. After 5000 seconds, the upstream deposition was 28 cm and the downstream deposition was 22 cm. Therefore, the average bed slope had increased from 1.23% to 2.09%. This can be one of the reasons why the average water depth in the numerical model was lower than in the physical model. The lower water depth also affected the average Froude number. Decreasing the water depth from 0.109 to 0.092 meters increased the average Froude number from 0.94 to 1.21.
The numerical model predicted the formation of laminae by small bed waves that formed on top of the antidunes. The initial bed waves emerged very small at the upstream end of the flume and grew in height and length as they moved in the downstream direction. The size of the bed waves affected the thickness of the layers they formed. At the upstream end of the flume, the height of the laminae was so small that they were not visibly resolved in the grid. In the middle of the flume, the layers were clearly visible in the figures and corresponded well with sizes observed in the field. However, further downstream, the bed waves increased in height and formed downstream-moving antidunes. Thereby, the sinusoidal bedding changed to a more complex sedimentary structure. Cross-bedding can be seen at the lower part of the strata. The water velocities were up to 1.8 m/s, and higher in the troughs than over the crests. This corresponded well with classical hydraulic theory and the Bernoulli equation. At 1000 seconds, the sixth wave from the upstream end had a strange shape, with a low velocity. This is an upstream-moving surface wave. All these observations will be discussed in more detail in the following section 'Upstream moving surface waves'.
Furthermore, Fig. 4 shows that the coarser particle sizes deposited at the upstream end of the flume.
This is a natural part of the sorting process in an alluvial channel, as the coarser material will have a lower sediment transport capacity. The sediments depositing at the downstream section of the flume therefore had a smaller average particle diameter and a higher sediment transport capacity. This is a likely reason why there was more deposition at the upstream part of the flume compared with the downstream part. Also, it explains why the bed slope increased. Note that the physical model study used a uniform sediment size, so this effect was only observed in the results from the numerical model. This is discussed in more detail in the following section 'upstream moving surface waves'. The upstream boundary condition is also a topic in the 'Discussion' section.

Sinusoidal structures from antidunes
A sedimentary structure from deposition under stationary sinusoidal antidunes was produced in the middle of the modelled geometry. Fig. 5 shows a more detailed view of this section at 5000 seconds, the end of the computation. The greyscale shows the d 50 of the grain sizes at the bed.

Bed waves
Looking carefully at the interface between the water and the bed in Fig. 5, it is possible to see some undulations on the bed, with a much smaller amplitude and wavelength than the antidunes.
These shapes are not reflected in similar forms at the water surface. The small waves move in the streamwise directions, as shown in Video 5 (Olsen, 2021). Bridge & Best (1997) observed similar processes in the laboratory and gave them the name low-relief bed waves. In the following, they will be called bed waves. They were also observed in the laboratory experiments by Núñez-González & Martín-Vide (2010) and Alexander et al. (2001). Other laboratory studies (Lunt & Bridge, 2007) also observed strata attributed to trains of small dunes forming on top of open framework gravel bars.
The numerical model gives more information about how the bed wave moved and formed the layers in the sedimentary structure. A detail of the flow field around a bed wave is given in Fig. 6.
The yellow arrow points to the front of the wave. The velocity vectors show a recirculation zone in front of the bed wave, with low velocities (white arrows). The return current stopped all bedload moving over the bed wave. The turbulent kinetic energy in Fig. 6 shows that the bed shear stress was very low in front of the wave. This caused sediments in front of the wave to deposit. At the top of the wave, the bed shear stress was relatively high, so that the finer sediments were eroded and coarser sediments dominated the bed grain size distribution. The bed wave therefore deposited a layer which had coarser sediments on top and finer at the bottom. This is also seen in the

Cross-stratification and downstream-dipping bedding
The lower part of the sedimentary structure in Fig. 5 shows downstream-dipping bedding at a high angle, forming cross-stratification. Since the angle of the bedding is different than the overlying deposits, it is likely that this structure was not formed by the same process as described previously. Fig. 8 shows more detailed longitudinal profiles of this structure during its formation. Video 8 (Olsen, 2021) shows the same process. The generation of the structures was caused by the growth of an antidune that moved in the downstream direction. There was a large recirculation zone behind the bedform, which stopped all the bedload. Therefore, the movement of the antidune must have been caused by sediment deposition at the top of the lee side, and then the forward movement had to be caused by sediment slides. This also explains why the structure has parallel lines equal to the angle of repose for the sediments. The movement of the front of the antidune followed a similar mode as dunes in rivers with low Froude numbers (Best, 2005). This was also observed by Spinewine et al. (2009): "Whereas upstream-migrating antidunes tend to be symmetrical, downstream-migrating antidunes tend to be asymmetrical, with a steep lee face similar to dunes". The steep lee side on a downstreammigrating bedform indicates movement by slides, which will produce downstream-dipping structures at an angle of repose for the sediments.
The formation of the structure with the downstream-dipping bedding took place from around 60-100 seconds after the start-up of the experiment. During this time period, the antidune migrated downstream. Afterwards, the antidune became stationary, and bed waves deposited sinusoidal laminae over the downstream-dipping bedding, resulting in cross-stratification. The equilibrium length of an antidune is related to the water depth (Kennedy, 1960). At time 67.5-90 seconds in Fig. 8, the length of the bedform was shorter than the equilibrium value. The recirculation zone at the stoss side was fairly large, and the sediments flowing over the crest deposited in slides. As the length of the antidune increased, the water flowing over the crest obtained a longer distance available for turning in the downwards direction. The recirculation zone then became shorter and smaller, as shown at 97.5 seconds in Fig. 8. Afterwards, the antidune was stationary in the horizontal direction and further deposition followed the wavy pattern given in Fig. 5. The downstream dipping bedding was formed as the antidune moved forwards, generating a cross-stratified structure with the overlying wavy deposits.
The thickness of a lamina is usually defined to be under 1 cm (Collinson & Mountney, 2019).
The horizontal size of the grid used in the present study was 1 cm. It was therefore difficult to compute lamina with a significant vertical slope, for example the downstream-dipping bedding in Fig. 8. In the present study they are called bedding, but with a finer grid it may be possible that the resulting size would be smaller than 1 cm. However, the resolution of the grid in the vertical direction is less than 1 mm, allowing the computation of thin near-horizontal structures that can be called laminae

Upstream moving surface waves
The numerical model had earlier been tested on downstream-moving antidunes (Olsen, 2016(Olsen, , 2017. However, also upstream-moving surface waves were produced by the numerical model in the current study. Cartigny et al. (2014) described a process called antidune breaking: "As the upstream flank of the antidune wave became oversteepened, flow over the antidune crests started to slide back against the incoming flow, producing rollers or breaking waves that migrated upstream as positive surges".
A snapshot of the surging surface wave is given at 1000 seconds of the 6 th antidune in Fig. 4. Fig. 9 and Video 9 (Olsen, 2021) show more details of a similar upstream-moving wave. Fig. 4 shows that these waves often occurred in the mid-section of the flume where the antidunes started to move downstream.
The irregularity of the bed at this location caused local variations in the water level, with accompanying changes in the pressure forces along the channel. An imbalance between upstream and downstream forces on a section of the flume could then cause a surge in the surface moving upstream (Cartigny et al., 2014). Fig. 9 shows the mid-section of the flume when the surface wave passed in the upstream direction. There are four bedforms shown in the figure. Two of the water waves moving upstream are indicated with yellow arrows and a letter. The downstream wave A moves first, followed by wave B.
This is seen directly on the water elevations in the figure, but also because the water velocity at the surface decreased as shown with the blue colour in Fig. 9. The velocities over the third bedform increased and there was some erosion at the bed, lowering the bed level. At 214 seconds, a second upstream-moving surface wave is seen over the second and third bedform. At 216 seconds, also the surface of the second antidune became skewed, and a region of low velocities emerged. Kennedy (1960) called the upstream movement of more than one surface wave for 'multiple breaking'. Momentum from the most downstream wave is transferred upstream, causing the antidune in Fig. 9 to move out of phase. Kennedy (1960) wrote ".. the breaking wave often rushed upstream onto the adjacent wave, thus precipitating its breaking". In the example shown in Fig. 9, the surge did not propagate fully to the most upstream (left) antidune.
The movement of the upstream wave was reduced a few seconds later, and at 228 seconds, the two upstream in-phase waves had regained their equilibrium shapes. However, both the surface and bed levels of the two most downstream antidunes had a lower amplitude than before. The upstream moving surface wave caused erosion at the crests of the bedforms, and left sediment deposits in the troughs. Alexander et al. (2001) observed that bed waves were important in the filling of the throughs: "Commonly, rapid erosion of an antidune crest leads to formation of a solitary, asymmetrical bedwave, up to a few centimetres high with a downstream face approaching the angle-of-repose. This bedwave migrates rapidly downstream into and fills, the trough area of the eroded antidune.". This is seen in Video 9 (Olsen, 2021). Fig. 9 shows four antidunes. The two most upstream bedforms (left) reached their equilibrium length of around 50 cm. The two most downstream antidunes were shorter than the equilibrium length at 200 seconds. During the upstream surface surge, the turbulence increased. This was partly due to the velocity gradients created by the surge and partly that the short bedform created an extra energy loss compared with the equilibrium sized antidune. The increased turbulence caused additional erosion over the third and fourth bedform. At the end of the surge, the erosion and deposition had extended the third antidune so that its length had reached the equilibrium value. The upstream-moving surface waves therefore contributed to a temporarily increase of the turbulence equilibrium in the train of bedforms.  The sedimentary structures in the downstream part of the flume at the end of the 5000 second simulation is shown in Fig. 10. As in the previous figures, the sedimentary structures were made up of strata with different grain size distributions. The thickness of the layers varied from 0.2 mm to 2 cm.

Sedimentary structures from downstream migrating antidunes
The longest structure was 30 cm. Due to the greater variation of the water depth in this downstream part of the flume, there was also a larger diversity in the thickness of the bedding/laminae.
The sedimentary structures did not have the same regular wavy pattern as in the mid-section of the flume.
The reason is that the antidunes at this location were not stationary, but migrated downstream.
Also, multiple occurrences of upstream-moving surface waves disrupted the formation of regular strata.
The steepest angles in Fig. 10 are between 10 and 12 degrees, both for the upstream and downstream dipping sets. This is lower than for the wavy structure shown in Fig. 5. The layers of the structure in Fig. 10 were formed by larger bed waves that eroded the sediments in front, as shown in Fig. 6 and Videos 5 and 8 (Olsen, 2021). This erosion removed all the structures except the layers deposited close to the troughs of the in-phase waves. Laminae formed at steep angles therefore did not remain.
The basic pattern-forming process producing the sedimentary structures in Fig. 10 is formation of convex-upward laminae with trough-shaped base deposition under antidunes. Two clear occurrences of convex-up layers are marked in the figure. As the antidune moved, bedwaves formed upstream and downstream dipping laminae. This is also indicated in Fig. 10. The occurrence of upstream-moving surface waves occasionally eroded the bed at the crest of the antidunes. This is a likely reason why there are no concave-upward layers in Fig. 10. The upstream-moving surface waves cause the troughs to be filled, which is indicated in the figure.
The water level and the bed in Fig. 10 are in phase, showing antidunes. Downstream dipping bedding is seen in front of each bedform. This would potentially create sets of cross-lamination. However, the deposits are quickly eroded by the downstream-moving antidune, so that the sets are fairly thin and it is difficult to observe any cross-lamination. The grid size of the bed is also 1 cm, preventing a high resolution of the downstream-dipping bedding and formation of possible sets of cross-strata.  Although Fig. 10 looks complex, the numerical model can explain every detail. An example is how the upstream dipping lamina ud1 formed. This is shown in Fig. 11 at three different time steps.
A medium-sized bed wave caught up with a downstream-moving antidune. The bed wave had a recirculation zone in front that caused erosion and could possibly have removed the crest of the antidune. However, for this particular bed wave, the turbulent kinetic energy was not sufficiently large to completely erode the bedform. This was partially due to the coarse material forming the crest at the bed of the antidune. The result was that the bed wave deposited a layer at an upstream dipping slope.
The set therefore dips upstream, although the antidune and the bed wave forming the structure moved downstream.

Discussion
The wavy sedimentary structure observed in Fig. 5 was caused by deposition under stationary antidunes. A question is how often this situation will occur in nature? If an antidune moves upstream or downstream, most of the sinusoidal laminae at the bed will be erased due to erosion (Kennedy, 1960). The result is a more complex structure like shown in Fig. 10. In the present study, the hydraulic parameters were chosen so that the Froude number was close to 1, a requirement for undular bedding to form. In nature, hydraulic conditions like channel width, water depth and channel slope will have a range of variation, both in time and space (Fielding, 2006).  The hydraulic parameters chosen for the present study replicate a laboratory study with downstreammigrating antidunes. However, the upstream boundary condition in the numerical model fixed the position of the most upstream antidunes, allowing sinusoidal laminae to form. Duller et al. (2008) observed irregular bedsets in a field study where a sedimentary structure formed after a glacial outburst flood. This was associated with a local hydraulic jump developing due to an isolated feature in the bed geometry. A natural river environment may include bed boulders or lateral contractions of the flow field causing similar conditions. It is therefore possible that isolated features could cause a wavy sedimentation pattern to develop also in nature. The sinusoidal laminae generated in the present study only formed for a few wavelengths.
The present study showed a breaking surface wave moving in the opposite direction of the flow ( Fig. 9). This has been reported from many laboratory studies (Kennedy, 1960;Alexander et al., 2001;Cartigny et al., 2014). However, it was not observed in the physical model of Núñez-González (2012), the study giving the input data for the present numerical model. There are several possible reasons for this.
The laboratory experiments of Núñez-González (2012) used a uniform sediment size, whereas the numerical model simulated a non-uniform sediment grain size distribution. Sediment sorting caused more sediment accumulation at the upstream end of the flume in the numerical model, leading to a steeper bed slope. The average Froude number therefore increased from 0.94 to 1.21. Kennedy (1960) observed that the breaking of the water surface over antidunes increased with the Froude number.
Another reason for the different upstream-moving surface waves could be numerical modelling errors.
Looking at Fig. 9, the upstream-moving surface wave has high turbulent kinetic energy close to the water surface. High turbulence levels at the water surface causes air entrainment, which was also seen in the laboratory study of Núñez-González (2012) and other physical model studies (Kennedy, 1960;Alexander et al., 2001;Cartigny et al., 2014). The surge of the upstream-moving wave was caused by momentum gradients in the channel. When air is entrained in the water, the density of the fluid will decrease.
The momentum of the surge is proportional to the fluid density and will therefore also be reduced.
The inability of the numerical model to include air entrainment may therefore cause the upstreammoving water surge to be larger in the numerical model than in the laboratory study. Technically, it is possible to include air entrainment in a CFD model (Almeland, et al., 2021). This is, however, a complex topic that is left for future studies.
A third reason for the overprediction of the upstream-moving surface wave is three-dimensional effects.
Bedforms observed in nature often have a more complex geometry than a two-dimensional model can reproduce (Fraccascia et al., 2016). Observations and videos from the flume study of Núñez-González (2012) showed multiple parallel surface and bed waves in the lateral direction, which were not always completely in phase. Turbulent diffusion in the lateral direction would cause a damping effect of the in-phase waves. Ideally, this would be best modelled with a 3D approach. Such computations would require substantially more computational time than what was required in the present study. This is also a topic for future research.
The average sediment size used in the present study was around 1.5 mm, which is several times larger than previous similar studies (Alexander et al., 2001;Cartigny et al., 2014;Vellinga et al., 2018).
This meant that the bed forms were less prone to erosion and therefore more stable. Ono et al.
(2020) used a base sediment of 0.14 mm and added some fraction of coarse material. Very complex irregular sedimentary structures emerged, which were not repetitive in the longitudinal direction.
However, similar structures were found in the field. The coarser material used in the present study could have caused less erosion during the upstream-moving surface waves. Kennedy (1960) reported that the bed levels were flat after passing of the wave. In the present study, the crests of the antidunes were eroded and the troughs were filled, but there were still some remains of the bedform left afterwards ( Fig. 9). This could be due to the relatively large size of the particles in the present study.
A common question in CFD analysis is the resolution of the grid. In the current study, the sediment bed was resolved with 1000 cells over a in the vertical thickness of 25 cm. Sedimentary structures extending mostly horizontally would then be well resolved. However, the grid size in the longitudinal direction was chosen to be 1 cm. When the sedimentary structures sloped with a large angle compared with the horizontal direction, the grid resolution will not resolve very large gradients. This is seen in Fig. 8, where the downstream dipping bedding shows a coarser stratification. A reduction of the grid size in the horizontal direction would significantly increase the already long computational times. It is therefore left for future studies. Wavy laminae were shown to be formed by smaller sediment waves moving over the bed of stationary antidunes (Fig. 6). A sorting process including variations in the bed shear stress and a small recirculation zone in front of the wave caused a coarse and a fine sediment layer to be formed simultaneously. A successive deposition from a number of these waves over a sinous-shaped antidune bed lead to an undulating sedimentary structure. The numerical model also showed how backsets can be formed in a similar process over downstream-moving antidunes (Fig. 11).

Conclusions
Downstream-dipping bedding at an angle of repose was found when a bedform with a significant vertical size moved downstream (Fig. 8). The numerical model predicted the formation of a recirculation zone at the lee side of the beform, together with a reduction of the bed shear stress in this location. The bedload was therefore stopped at the top of the lee side of the crest, causing sediments to deposit and slide down on the lee side. The numerical model was able to compute the resulting bedding from this process, although with a coarse resolution.
The numerical model includes equations for predicting the turbulent kinetic energy of the flow, which can be used to predict the bed shear stress. Video 9 (Olsen, 2021) shows how antidunes with equilibrium dimensions have much less turbulent kinetic energy than bedforms that are shorter.
An antidune not in equilibrium will therefore have a higher bed shear stress and more potential to change its dimensions. The CFD model thereby gives an explanation why equilibrium dimensions of antidunes exist, and thereby also why their geometrical parameters can be predicted by empirical formulas.