Calculation of dune morphology - Wiley Online Library

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 110, F04013, doi:10.1029/2004JF000171, 2005

Calculation of dune morphology Søren Tjerry DHI-Water and Environment, Portland, Oregon, USA

Jørgen Fredsøe Technical University of Denmark, Lyngby, Denmark Received 27 May 2004; revised 14 June 2005; accepted 19 August 2005; published 30 November 2005.

[1] The shape and dimensions of dunes formed in an erodible bed and exposed to flowing

water are investigated by using a numerical model. The flow is calculated using a curvilinear flow model, where we combine a k-e turbulence closure scheme with a sediment-transport description and the equation of continuity for the sand bed. The complete flow description allows for a discussion about what determines the length of dunes: From the equation of continuity for sediment, it turns out that the sediment transport over dunes is proportional to the local height of the dune. This means that the maximum height of the sand wave occurs at the location of maximum sediment transport. We relate the dune length to the location of maximum sediment transport. It is demonstrated that the effect of streamline curvature is essential for determining the location of maximum sediment transport at low bed shear stresses, due to an overshoot in bed shear stress. For large dimensionless bed shear stresses the downstream decay of turbulence from the former dune is demonstrated to play a major role. Citation: Tjerry, S., and J. Fredsøe (2005), Calculation of dune morphology, J. Geophys. Res., 110, F04013, doi:10.1029/2004JF000171.

1. Introduction [2] A plane streambed is usually unstable and tends to break up with different kinds of sand waves formed [e.g., Guy et al., 1966]. One of the most important configurations of sand waves are the dunes, which are formed in the subcritical regime and sometimes related to the presence of a hydraulically rough bed [Engelund and Hansen, 1973]. While the formation of sand waves in alluvial streams has attracted much attention over several decades [e.g., Smith, 1970; Kennedy, 1971; Reynolds, 1976; Engelund and Fredsøe, 1982; McLean, 1990; Fredsøe, 1996], very little research has focused on modeling the fully developed sand waves. This is despite its significance to engineering, since dunes have a great impact on the flow resistance and the sediment transport. Empirical work on dune properties has been done among others by Van Rijn [1984], who suggested empirical relations among height, length and steepness of dunes. Since these relations are widely used, comparison with Van Rijn’s empirical curves is included in this paper. No complete theoretical understanding of the properties of dunes is available. [3] A number of previous studies have formulated theoretical models for the flow over dunes. McLean and Smith [1986] used the law of the wake to describe the flow past fully developed bed forms, including the effect of flow separation from the upstream dune. Nelson and Smith [1989] refined the model by including the total influence Copyright 2005 by the American Geophysical Union. 0148-0227/05/2004JF000171

from all upstream dunes. Mendoza and Shen [1990] used an algebraic Reynolds stress closure to describe the flow over bed forms using a Cartesian grid. Calculations with curvilinear models have been performed by, for example, Van der Knaap et al. [1991], Richards and Taylor [1981], Zijlema et al [1995] and Yoon and Patel [1996]. Very few scientists have made the additional step to link the flow together with a sediment-transport description to calculate the shape of self-formed sand waves. One of the first was Exner [1925], who used a depth-integrated flow description to predict the asymmetric evolution of the dune shape, but the model was unable to describe the equilibrium situation, where the height becomes constant in time. Smith [1970] used a constant eddy viscosity concept to calculate the flow over sinusoidal bed forms, and included a discussion of the properties of fully developed bed forms. The measured flow behind a sharp discontinuity in the bed geometry was used to discuss the flow properties further downstream on the plane bed. Smith observed a peak in the bed shear stress a certain distance downstream of the discontinuity, and related the length of the bed forms to this local peak. Fredsøe [1982] also used the flow behind a backward facing step to discuss the properties of sand waves. A full morphological calculation of the dune shape was made based on a semi-empirical flow approach, where Fredsøe combined a depth-integrated flow description together with a description of the flow relaxation behind a step. This was based on the measurements by Bradshaw and Wong [1972], who described the recovery of the flow after the disturbance due to the step as a relaxation process. Fredsøe demonstrated that a detailed description

F04013

1 of 13

F04013

TJERRY AND FREDSØE: CALCULATION OF DUNE MORPHOLOGY

F04013

H is the dune height as shown in Figure 1a. An expression for the variation of the sediment-transport q along the dune can be obtained from the sediment transport equation (also called the Exner equation), @q @h ¼ : @x @t

ð3Þ

From equation (1) we get a

Figure 1. (a) Definition sketch and (b) migration of dunes. of the sediment-transport pattern was essential for the calculation of dune properties. At low dimensionless bed shear stresses, the effect of the local slope on the bed load transport was shown to have an important impact in predicting the dune length, while at high dimensionless bed shear stresses, the impact of suspended sediment had a similar important effect for the prediction of the length. [4] The work by Fredsøe [1982] (hereinafter referred to as JF-82) is taken as the starting point for this paper. The semi-empirical flow description is replaced by an accurate numerical flow solution, based on turbulence modeling, mainly using the common k-e model. One shortcoming of the JF-82 model was that the effect of flow curvature was not included, meaning that the JF-82 model only can describe a ‘‘train’’ of sand waves, since the shape is calculated by looking at the relaxation of the boundary layer occurring downstream of the former dune crest. In nature, solitary sand dunes are also observed, [e.g., Termes, 1984; Engelund, 1971], which suggests that relaxation is not the only mechanism to explain the length of sand waves. The present work proposes that the effect of streamline curvature is an additional controlling mechanism for the geometry of fully developed sand waves. The effect of curvature has previously been included in the linear analyses of flow over a sinusoidal bed [e.g., Benjamin, 1959; Smith, 1970; Engelund, 1970; Fredsøe, 1974; Richards, 1980]. Engelund [1970] and Smith [1970] applied this in combination with an eddy-viscosity model to determine the preferred wavelength in the initial phase of growth from an initial plane bed. However, the physics in the flow is, although related, very different in the linear case and in the present strongly nonlinear case.

@q @h ¼a : @x @x

q ¼ q0 þ ah;

a¼

qD : H

ð2Þ

ð6Þ

where q0 is the sediment transport in the trough of the dune where h is zero (compare Figure 1). Equation (6) shows that the local sediment transport is proportional to the local height h of the dune.

3. Governing Equation for the Bed Form Shape [6] The morphological calculation in the model is as follows: It is assumed that the sediment transport occurs as bed load only; that is, it responds immediately to changes in the spatial shear. The inertia of the bed load particles is ignored. Assuming bed load to be the only sediment transport mode, the constant q0 appearing in equation (6) will be zero because no sediment transport takes place in the trough where h is equal to zero (see Figure 1). Now equations (2) and (6) give h q ; ¼ H qtop

ð7Þ

where qtop denotes the sediment transport on the dune crest. The bed load transport can be calculated by the Meyer-Peter and Mu¨ller [1948] formula, which reads ð8Þ

where q is the Shields parameter defined as q¼

ð1Þ

The migration velocity is determined by the amount of deposited sediment qD (volumetric, including porosity n) that is trapped on the rear-side of the dune (see Figure 1b). The migration velocity is determined by

ð5Þ

By integration of this equation, it is found that

q 8ðq qcr Þ3=2 F ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ; ð1 nÞ ðs 1Þgd 3

[5] We consider 2D uniform sand waves migrating with the migration velocity a in the downstream x-direction without changing its shape h with time t, i.e.,

ð4Þ

This can be inserted into equation (3) to give

2. Variation in Sediment Transport Along the Dune Surface

hð x; t Þ ¼ hð x atÞ:

@ @ ¼ : @x @t

U*2 ; ðs 1Þgd

ð9Þ

where U* is the friction velocity, s relative sediment density, d mean grain diameter, g acceleration of gravity and qcr the critical Shields parameter above which the sediment grains begin to move. Meyer-Peter and Mu¨ller [1948] take the critical Shields parameter to be 0.047, which also is being used in this paper neglecting the slight grain size dependency of qcr. On a bed with a longitudinal slope, q in equation (8) must be modified to q*, which includes the

2 of 13


F04013

F04013

empirical way. An empirical expression for the bed shear stress variation due to the flow separation was included based on bed shear stress measurements on a horizontal bed downstream a backward facing step with a height Hs, mainly based on the measurements by Bradshaw and Wong [1972]. This expression is given by cf x x ¼ 1 exp 0:377 5 1:06 ; cf 0 Hs 350Hs

ð12Þ

where cf is the skin friction coefficient and cf0 is the maximum value of this coefficient, which occurs around 16 step heights downstream of the step. The behavior of cf according to equation (12) is shown in Figure 2. The cause of the local maximum in cf is not totally clear. Other flow studies of the same problem do not as clearly reveal this local maximum, which might depend on properties like bed roughness and Reynolds number. Without this local maximum, the JF-82 model cannot predict a maximum in the bed shear stress over a dune. The variation in the bed shear stress tb is related to cf by tb ¼ rcf V 2 :

Figure 2. Flow downstream of a backward facing step. (top) Picture of streamlines calculated using the RSM model. To illustrate the separation zone, the streamline interval is smallest at the top of the step. (bottom) Comparison of the predicted variation in the skin friction coefficient with the experiments by Bradshaw and Wong [1972] and the semi-empirical expression by Fredsøe [1982]. forces acting on the bed load particles due to the longitudinal bed slope [see, e.g., Fredsøe and Deigaard, 1992, p. 205] q* ¼ q m

@h ; @x

ð10Þ

where m is a dimensionless constant equal to qc/tan (fs). Here fs is the angle of repose for sediment, giving m a value around 0.1. The factor (1-n) appearing in equation (8) is introduced as the volumetric transport rate including pores is applied in this paper. Inserting equation (8) into equation (7), the following first-order differential equation in determining the shape h is obtained: dh ¼ dx

2=3 0:047 qtop h q 0:047 þ : m H m

ð11Þ

Equation (11) is the governing equation, which connects the shape h to the spatial variation in bed shear stress.

4. Flow Description Along the Dune [7] The flow over dunes is highly complex due to its non-uniformity in depth and the flow separation. In the earlier JF-82 model, both features were included in a semi-

ð13Þ

Here r is the fluid density and V is the depth-integrated average flow velocity, which varies along the dune due to changes in the bed level and thus the water depth. This variation is included applying the continuity equation for the fluid. For a plane water surface occurring at low Froude numbers, this will modify V in the following way: V ð1 h=DÞ ¼ Vtop ð1 H=DÞ:

ð14Þ

Vtop is the velocity at the dune top, where h = H, and D is the mean water depth along the dune (compare Figure 1a). The variation in cf due to equations (13) and (14) is called the cross-section effect, since the variation is due to differences in the cross-sectional area. This suggests that the Shields parameter can be rewritten in terms of the relaxation function applying equations (13) and (14). q ¼ qmax

cf ð xÞ 1 H=D 2 ; 1 h=D cf 0

ð15Þ

where qmax is an input value for the Shields parameter in the calculations. [8] The simple flow model described above does not include any effect of the bed (and hence flow) curvature on the bed shear stress distribution. The convex curvature along the stoss part of the dune gives additional acceleration of the fluid resulting in an increase in the related bed shear stress. This increase is quite similar to the local increase in the flow velocity in the well-known flow case around a circular cylinder. The increase is largest at that location where the bed curvature has its maximum, so the effect decreases further downstream of the dune as sketched in Figure 1. Owing to the contraction of the flow along the gentle part of a real dune, a faster recovery

3 of 13

F04013


F04013

chosen to include only a very short description of the model in Appendix A.

5. Verification of the Hydrodynamic Model 5.1. The Backward Facing Step [9] The flow over dunes changes from flow exposed to a strong adverse pressure gradient just downstream of the dune crest to a slightly favorable pressure gradient along the gentle positive slope of the dune. Downstream of the dune crest, the flow separates and the slope of the sand wave is here nearly equal to the angle of repose, indicating that the flow in the separated zone is quite weak. This part of the flow has some similarities to the flow behind a backward facing step, which is one important feature which the flow model must be able to describe. Figure 2 shows the result of the calculations of the flow past a backward facing step. In the upper part of the figure, the streamline pattern is shown. The bottom part of Figure 2 presents a comparison between the different model predictions of the bed shear stress tw (through the variation in cf) which in the turbulence models are calculated from the non-equilibrium condition at the wall, u ¼

Figure 3. Predicted time-averaged flow velocity: dashed: RSM, solid: k-e model. (a) Horizontal flow velocity. Measurements (circles) by Nelson et al [1993]. (b) Like Figure 3a but with experiments by Bennett and Best [1995]. (c) Vertical flow velocity, with measurements by Nelson et al. [1993].

of the boundary layer must be anticipated in this case as compared to the case of flow over a horizontal bed downstream backward facing step. In order to obtain a complete flow description, including all the above mentioned flow features, the flow in this paper is calculated by computational means, using a standard turbulence model. Because of the flow separation, a two-equation k-e model was used rather than a one-equation k-model (see Launder and Spalding [1972] for a review of turbulence models). Also, the more sophisticated differential Reynolds Stress Model (RSM) was applied in a few cases to gain insight into some specific physical matters and to demonstrate the sensitivity in the results due to the choice of sophistication of the turbulence closure scheme. Since these models today are standard tools [see, e.g., Wilcox, 1993], we have

tw 30y ; pffiffiffiffiffi log k tk kN

ð16Þ

where k is the Von Karman’s coefficient, y the distance from the bed, kN the bed roughness and tk is the shear stress due to turbulent kinetic energy k. k is related to pffiffiffiffiffi the shear velocity by k = U*2/ cm , where cm is the return-to-isotropy constant (0.09). Figure 2 shows that the k-e model is incapable in reproducing a local maximum, while the differential RSM applied in the present case actually produces such a maximum. The possible presence of a local shear stress maximum can be related to the structural changes of the turbulence described by the RSM model. Also included in Figure 2 is the effect of a downstream triangular contraction used to simulate the presence of the downstream dune. The separation zone becomes shorter in this case due to the contraction. 5.2. The Flow Over a Dune [10] The experimental works by Nelson et al. [1993] and Bennett and Best [1995] are used for verification of the numerical models. A later experimental work by Venditti and Bennett [2000] could also be used, but it has fewer measured sections along the dune. Nelson et al. [1993] used a LDA to measure the velocity and turbulence in the flow over concrete dunes shaped by a cosine-shaped stoss and a linear lee side. The water depth was 21.5 cm and the experiment considered a nondimensional dune length L/D = 3.36 and dune height H/D = 0.168. The hydraulic roughness kN was not reported, and the relative roughness kN/D is estimated to be 0.002 (J. Nelson, personal communication, 2002). Bennett and Best [1995] covered their dune with glass spheres with a diameter of 0.22 mm, the roughness estimated herein to be twice this diameter. The water depth was 14 cm, and the dimensionless dune length and height were 4.5 and 0.286 respectively, so the dune considered by Bennett and Best is higher and longer

4 of 13

F04013


Figure 4. Decay of turbulence intensity k downstream a rearward facing step and along a dune. than that considered by Nelson et al. In the calculations presented in Figure 3, the upstream dune crest is represented by a backward facing step. Figures 3a and 3b show the comparison between measured and predicted horizontal velocity fields. Run 4 from the Nelson et al. experiments was considered, where seven verticals along a dune were measured. The general agreement is satisfactory, but there is a tendency toward an under prediction of the horizontal velocities at the bed. This can be partly tuned by changing (decreasing) the bed roughness. Figure 3c shows a similar comparison of the predicted vertical velocities with the Bennett et al. measurements and satisfactory agreement is also found here. The differences between model predictions of the flow over dunes and downstream a backward facing step is illustrated in Figure 4 by comparing the decay of the turbulence in the water column. The measurements by Mueller and Robertson [1963] are considered for the backward facing step, while the Nelson et al. [1993] data are considered for the flow over a real dune. In the case of the backward facing step, the tendency is that the models yield a faster decay in the turbulence. The flow past a real dune shape does not encounter the same problem because of the contraction of the flow lines and the associated faster decay of the turbulence.

F04013

shape of dunes calculated by using a k-e model and a RSM model are shown in Figure 5. In Figure 5a, the Shields parameter is only slightly above the critical value (0.047 in the Meyer-Peter formulation) and both turbulence models predict very long bed forms: the k-e model predicts a sand wave steepness (height to length) of approximately 1:70, while the similar result from the RSM is 1:90. In Figure 5b the Shields parameter has been increased to be more than twice the critical value, and the steepness is increased significantly, namely to 1:25 and 1:33 respectively. In all of the calculation cases, the dune stops when the bed level begins to decline, as the bed shear stress in a divergent flow will drop immediately with flow separation as a result. A relatively sharp corner is formed at this location on the dune, and this requires a numerical grid with large curvature. Such a grid causes very high bed shear stresses locally as discussed above in section 4, and also experienced by Van der Knaap et al. [1991]. For this reason, the downstream steep dune front is not included in the numerical scheme for the flow calculation. Instead, the flow is calculated over the predicted dune until the crest has been reached. Downstream of the crest, the dune is continued with a straight line at the same elevation of the crest level. By avoiding the downstream front in the numerical solution, the bed shear stress is still satisfactorily described since the presence of the separation zone on the lee side of the dune forces the flow outside the separation roller to continue more or less tangential to the slope of the sand wave at the crest, which can be seen in the detailed photos from Guy et al. [1966]. The results presented in Figure 5 are in line with the earlier findings in JF-82, who also predicted that the steepness would increase with increasing q at low q-values. The interesting thing is that it now can be explained differently, namely as a ‘‘curvature-effect’’: As shown in the depicted cf variation in Figure 5, the curvature of the dune-surface causes a local maximum in the bed shear stress which, unlike the plane bed case (Figure 2, bottom), is even predicted by the k-e model. This maximum is therefore not entirely determined by the relaxation process as suggested in JF-82. Another contribution to the presence of a maximum in the bed shear

6. Calculations of the Bed Form Shape [11] As a first approximation, the flow past a backward facing step is solved to calculate the initial value of cf (x) in an iterative procedure. Next the shape is calculated from equation (11), and a new cf (x) is calculated from the hydrodynamic code with the new shape inserted as the bottom boundary, where the curvilinear grid is updated to fit the modified boundary. Equation (11) is solved by a simple Euler scheme and the calculated shape is filtered after each iteration applying a numerical filter [see, e.g., Richards and Taylor, 1981; Jensen et al., 1999]. The standard filtering operator used here is the same as equation (12) of Jensen et al. [1999]. Without using such a filter, the solutions become unstable and will diverge. Appendix 1 of Jensen et al. [1999] describes the whole morphological scheme applied in the present paper. Examples of the

Figure 5. Steady dune profile obtained by the k-e and the RSM model for (a) H/D = 0.053, qtop = 0.06 and (b) H/D = 0.15, qtop = 0.106. The skin friction variation is that predicted by the k-e model.

5 of 13

F04013


F04013

further purely a function of the Shields parameter, and therefore @q @q @q ¼ : @x @q @x

ð18Þ

The variation in q with respect to x is due to two effects as seen from equation (10): the ‘‘cross-section effect’’ and the relaxation. The latter seems to be insignificant based on tests using the computational code. Figure 6a shows an example on the calculated shear stress distribution over a dune with triangular shape. It is observed that the whole shear stress variation at the dune crest is due to the ‘‘cross-section effect,’’ while the skin friction coefficient has approached a constant at a distance of just 10– 12 times the dune height downstream the former dune crest. Figure 6b shows the similar picture for a more realistic shape, calculated by the JF-82 method. By neglecting the variation in the skin friction coefficient at the dune crest, we get by differentiation of equation (15) the following variation in the dimensionless bed shear stress at the crest where h is equal H: @q 1=D @h ¼ qmax : @x 1 H=D @x

Figure 6. Contributions to the bed shear stress past (a) a constant slope bed form and (b) a realistic dune shape calculated by the JF-82 method. The definition of the ordinate axis is shown on the individual curves.

By inserting this equation into equations (17) and (18) we get H D

1 HD

stress can instead be attributed to the fact that along a concave curved boundary, the negative (or favorable) pressure gradient increases with curvature, so a maximum occurs close to the dune crest, where the bed shape still is slightly convex. In combination with the high flow velocities near the crest, where the cross section (i.e., local water depth) is small, this effect contributes to the local maximum in the bed shear stress.

ð19Þ

¼

FD : 2q dF dq

ð20Þ

In Equation (20), we have replaced q and qD by their dimensionless versions F and FD. In the bed load case, it can be assumed that all bed load transported over the top will deposit on the front. Inserting the bed load formula equation (8) into equation (20) yields H=D q qc : ¼ 3q 1 H=D

ð21Þ

This expression is shown with the solid line in Figure 7.

7. Dune Height, Steepness and Length 7.1. Dune Height [12] As demonstrated above, a solution for the shape of a dune can be obtained for any specified value of H, and as the length of the dunes depends on their height, the prediction of this height is first discussed. Fredsøe [1979, 1982] and Deigaard and Fredsøe [1987] suggested a way to calculate the dune height by considering a moving front on which the trapped sediment in the separated region causes the dune to migrate forward. This method is reviewed in the following, and is discussed with the improved flow description available. In the case of pure bed load, the method uses the continuity equation (5) that can be rewritten into @q qD @h ¼ ; @x H @x

ð17Þ

where qD is the deposited (trapped) sediment on the steep stoss part of the dune. At the flat, nearly horizontal part of the dune near the top, the sediment transport is

Figure 7. Dune height predicted by equation (21). The data are from Guy et al. [1966].

6 of 13


F04013

F04013

periodic with a wavelength equal the dune length. This last condition is obtained by the following iterative procedure. First, at the inlet of the computational domain, a Vanoni profile is applied for the vertical distribution. Further development is calculated from the transport equation (23). The profile at the outlet is inserted at the inlet as a new input profile, and the calculations are repeated until periodicity is obtained. The amount of trapped sediment shown in Figure 8a is found as that fraction, a, of the suspended sediment that enters the separation zone downstream the crest. This depends on the Rouse parameter z, defined by z¼

Figure 8. Behavior a backward facing (b) Calculated lines for the case H/D = relative variation in next is constant.

of suspended sediment behind step. (a) Sketch of deposition. of constant sediment concentration 0.2, D/d = 2500 and z = 2. The concentration from one line to the

[13] At higher Shields numbers, a large fraction of the sediment transport occurs as suspended sediment and in this case only a certain fraction a of the sediment transport over the dune crest is trapped on the stoss part. FD ¼ aFtop :

ð22Þ

kUf : ws

ð25Þ

Figure 8b shows an example of the calculated distribution of suspended sediment past a step. Very few measurements of the distribution of suspended sediment are available to confirm the predicted distribution. Venditti and Bennett [2000] did measure this distribution, but over a starving bed, so that spatial changes close to the bed due to deposition and pickup are not correctly reproduced. Our prediction of the concentration profile suggests a faster decrease in concentration than measured by Venditti and Bennett [2000]. They attributed this to macro-turbulent events, which will require a Large Eddy Simulation or Direct Numerical Simulation to include in the description. Such a detailed 3D flow model is at present too computationally intensive for detailed morphodynamic modeling. It is correct that a time-averaged model like

With the present detailed flow description, it is possible to estimate the value of a, which was taken to be equal one in JF-82. Further, the lag between the local dune height and the local transport of suspended sediment can be described in more detail with the present flow description. This lag is present as it takes some time for the suspended sediment to settle after being picked up from the bed. In the stability analysis, this lag is able to explain the transition regime from dunes to plane bed [Engelund and Fredsøe, 1982]. In the present analysis, the inclusion of the detailed description of suspended sediment is used to describe the slight decrease in dune height and the increase in wavelength at large Shields numbers. The transport of suspended sediment is described by assuming that the distribution is determined by advection, diffusion and settling as done by Engelund [1970]. The transport equation reads @c @ ðucÞ @ ðvcÞ @ ðws cÞ @ @c þ þ ¼ þ nt ; @t @x @y @y @y @y

ð23Þ

where c is the volumetric concentration of suspended matter, nt the eddy viscosity and ws is the settling velocity of the sediment grains. As boundary condition we use the bed concentration concept, cb ¼ cb ðqÞ;

ð24Þ

in which the formulation by Engelund and Fredsøe [1976] for cb is used. A zero flux boundary condition is used at the water surface, and finally the concentration profiles must be

Figure 9. Ratio between measured and calculated dune migration velocity. (a) Only bed load assumed trapped. (b) All sediment assumed trapped. (c) Model prediction. The data are from Guy et al. [1966].

7 of 13

F04013


Figure 10. Impact of selected dune height on predicted dune steepness. the k-e model gives a slight underestimate of turbulent kinetic energy and thereby the diffusion coefficient far away from the bed (see Appendix A). Close to the wall, however, k is quite well reproduced. [14] Spatial changes in suspended load due to sedimentation and pickup occur mainly in the near wall layer, while the transport farther away is nearly constant along the dune. Since our calculations are based on gradients in the sediment transport, the diffusion concept should therefore be appropriate for the present investigation. An indirect proof can be obtained by comparing the calculated prediction of trapped sediment on the dune front. This comparison is made from comparing calculated and predicted migration velocity of the bed forms as shown in Figures 9a – 9c. In Figure 9a, only the bed load is assumed trapped. In Figure 9b, all the sediment (bed load and suspended load) is trapped, while Figure 9c is based on the model calculations. In JF-82, only the bed load was assumed trapped, but it seems from Figure 9c that the migration velocity is slightly better described by adding a certain fraction of the suspended sediment; however, the scatter is large. Equation (19) is used to calculate the dune height as shown in Figure 7 in which F now is the sum Fb + aFs. This calculated result is compared with Guy et al.’s

Figure 11. Predicted variation in dune length for the bed load case. The data are from the 0.93-mm case of Guy et al. [1966].

F04013

Figure 12. Predicted position of maximum skin friction on a low dune. [1966] data, and even though the scatter is large, the prediction matches qualitatively with the experimental data. The predictions shown in Figure 7 agree with the observed trend that the dune height becomes smaller at high Shields numbers for fine sediment (large values of D/d). In the field a similar behavior is more difficult to observe which often is related to the large timescale for large dunes to change properties. One example of a decrease in the dune height due to an increase in the amount of suspended sediment is the Missouri River, where U.S. Army Corps of Engineers [1969] reported a drop in dune height and an increase in length due to a decrease in temperature and an associated increase in the amount of suspended sediment. 7.2. Length of Dunes [15] The predicted dune length is obtained by dividing the dune height with the steepness. The predicted steepness of the dunes is found by first calculating shape of the dunes as described above, and then relating the length to the location of the bed load transport maximum. Figure 10 shows the predicted variation in steepness H/L with q, which is the only dimensionless parameter appearing in the steepnessvariation at low Shields numbers. The predictions fit the coarse sand (i.e. no ripples) laboratory measurements by Guy et al. [1966] well. The predicted length and the steepness depend on the correct choice of the dune height. This is illustrated by Figure 10, where in one case the height has been calculated using the equation (20), and in the other case, it has been made equal to a constant, namely H/D = 0.1. Figure 7 reveals that for the Shields number larger than 0.07, the predicted dimensionless height becomes larger than 0.1. Figure 10 shows that if the dune height is assumed to be smaller than the predicted one the predicted steepness becomes overestimated (for q larger than 0.07). [16] For steep dunes, the separation-reattachment process is important for determining the wavelength because it is setting a lower value for this property. For lower dunes, the length seems more related to the curved shape of the dunes, and the dune length seems to be independent of the relaxation process at low Shields parameters, where the

8 of 13

F04013


Figure 13. Calculated dune shape in which suspended sediment is included in the model for (a) qtop = 0.4 and H/D = 0.4 and (b) qtop = 1.0 and H/D = 0.05. The definition of the ordinate axis is shown on the individual curves. height of the dune becomes very small (compare Figure 7). Figure 11 shows the prediction of the dune length. The present predictions suggest a quite constant dune length at low Shields numbers in agreement with the measurements. This improves the earlier findings JF-82, who reported a significant increase in length at small Shields numbers. This illustrates the importance of including the curvature effect in the domain, since the predicted length thereby becomes much longer. Also shown in Figure 11 is the empirically based suggestion by Van Rijn [1984] that the length L constantly is equal to 7.3 times the water depth. Figure 12 shows that the calculated position Lm of the maximum bed shear stress at small Shields parameters (or low dunes) occurs far downstream of the maximum that occurs downstream of a backward facing step. [17] At higher Shields numbers, the suspended sediment must be included in the morphodynamic description. Also,

Figure 14. Modification in the dune shape with increasing Rouse number z, corresponding to an increase in the amount of suspended load.

F04013

Figure 15. Dune steepness prediction including suspended sediment and comparison with field and laboratory data. the dune shape is modified owing to the presence of suspended sediment. Two examples of the shape calculations at low and high Shields numbers are shown in Figures 13a and 13b, respectively. It is shown that the suspended sediments cause a large downstream shift in the position of the maximum sediment transport, so the predicted dune becomes much longer. Figure 14 illustrates how the suspended sediment modifies the shape. Three calculation examples are shown with increasing values of the Rouse parameter z, and keeping the remaining parameters the same. The amount of suspended sediment increases quickly with z, so the moderate variation in z shown in Figure 14 corresponds to an increase in qs of more than 50%. For increasing amounts of suspended sediment (due to an increase in z), the shape becomes relatively steeper during its growth, followed by a very mild increase toward the dune crest. It is difficult to compare this with measurements, but Figure 14 is in agreement with the findings by Klaassen et al. [1986], who reported that the shape of dunes at high velocities was characterized by a fast increase in the bed elevation in the relaxation region, and then a slowly growing bed elevation that may continue for a very long distance. [18] Figure 15 shows the predicted variation in dune steepness with the Shields parameter for different values of the depth D. The dune height has been taken to be constant, and equal to 20% of the water depth to avoid too many variables in the input parameters. If instead, the predicted height from equation (20) is applied, any deviation from the curves depicted in Figure 15 can hardly be detected. The data are laboratory measurements from Guy et al. [1966] or field measurements partly provided personally

9 of 13

F04013


F04013

[21] The sediment transport calculations are in the present paper based on the time-averaged bed shear stress. Nelson et al. [1995] emphasized the importance of including a more detailed sediment transport description based on the local characteristics of the fluctuating part of the bed shear stress as done by Sumer et al. [2003]. This will locally lead to an increase in the sediment-transport rate and make the peak in sediment-transport more pronounced. This will certainly modify the sediment transport description of some dune heights downstream of the point of separation, and shall be incorporated in more detailed modeling of the dune behavior. Figure 16. Dune length prediction including suspended sediment. Data are from Guy et al. [1966].

by C. F. Nordin or described by Julien and Klaassen [1995]. It is seen that the predicted steepness lies centrally within the very scattered measurements. Also, the tendency that the fine sand will lead to less steep dunes corresponds to the measured trend. [19] The predicted dune length, as in the bed load case, was obtained by dividing the dune height (Figure 7) with the steepness (Figure 15). The prediction is shown in Figure 16. Also shown in the figure is the suggestion by Van Rijn [1984] that the length L is about 7.3 times the water depth. Consistent with the JF-82 theory and the data, the present theory predicts a certain increase in length with increasing Shields parameter.

Appendix A: Short Description of the Turbulence Models [22] The k-e model applies the eddy viscosity concept for describing the Reynolds stress, i.e., ui uj ¼ nt

@ui @uj 2 þ dij k; @xj @xi 3

where ui uj is the Reynolds stress tensor, xi is the coordinate, nt the eddy viscosity, ui the velocity vector, dij the Kronecker delta and k the turbulent kinetic energy. The last term in equation (A1) ensures consistency with the definition of the turbulent kinetic energy, and it is included by implicitly modifying the pressure. The eddy viscosity is given by nt ¼ cm

8. Conclusion and Final Remarks [20] The shape and dimensions of dunes have been treated theoretically using a turbulence closure to describe the flow. The basic morphological approach is similar to that originally developed by Fredsøe [1982], but the semiempirical approach is now replaced by an accurate flow description based on the model. For low Shields parameters, it is concluded that the length of dunes is determined by a combination of a relaxation process and a ‘‘curvatureeffect,’’ which causes a peak in the boundary shear over a real sand wave. The new approach mostly differs quantitatively from the original approach of Fredsøe at very low Shields parameters. The present approach yields a fairly constant dune length, which is in agreement with the observations, while the method of Fredsøe yields a dune length that approaches zero at the critical bed shear stress. The present model suggests that at very low Shields parameters the dune length become determined purely by the local curvature. Applying a transport equation for the suspended sediment includes the effect of the lag in the suspended sediment. The amount of trapped sediment on the steep lee side of the sand wave is calculated and shown to be important for the determination of the dune height as well as the dune length. It is shown that an assumption of local equilibrium for the suspended load leads to a substantial overestimation of the decrease in dune height at high Shields parameters. A modified calculation method for the dune height is outlined in which the influence of the suspended load on the dune height is significantly weaker than in the original approach of Fredsøe.

ðA1Þ

k2 ; e

ðA2Þ

where e is the dissipation of energy. The turbulent kinetic energy is found by solving the transport equation, ui

@k @ nt @k @ui @uj @ui ¼ þ e; þ nt @xj @xi @xj @xi @xi sk @xi

ðA3Þ

where sk is a Schmidt number for the turbulent kinetic energy. In the present work the turbulent kinetic energy is solved with inflow conditions at inflow boundaries, and zero flux everywhere else. A non-equilibrium law of the wall for the velocity profile accompanies the zero flux condition at the bed given by equation (16). The transport equation for the dissipation is given by ui

@e @ nt @e @ui @uj @ui e2 þ ce1 cm k ¼ þ ce2 ; ðA4Þ @xj @xi @xj k @xi @xi se @xi

where se is a Schmidt number for the dissipation, and ce1 and ce2 model constants. The boundary condition approach is the same as for the turbulent kinetic energy, except at walls where the dissipation is calculated by assuming equilibrium between production and dissipation of turbulent kinetic energy. The following parameters are used in the k-e model [Rodi, 1993]: cm = 0.09, sk = 1, se = 1.3, ce1 = 1.44, ce2 = 1.92. [23] In the Reynolds Stress Model (RSM) the three separate normal stresses are handled instead of just the turbulent kinetic energy, and the Reynolds stress is solved

10 of 13

F04013


as well. In the RSM there are thus five turbulence variables instead of the two in the k-e model. The principles are the same as for the k-e model. The Reynolds stresses are found from transport equations where the boundary conditions are the same as used for the k-e model, i.e., zero flux for the stresses, except at inlet boundaries, while the Reynolds shear stress is used in the non-equilibrium wall function for the flow velocity. For the dissipation, the equation is the same as in the k-e model, except that the production of turbulent kinetic energy is given by its exact expression. The applied RSM is the so-called Gibson-Launder version; see Launder et al. [1975] for details. [24] The turbulent kinetic energy is shown in Figures A1a and A1b. In analyzing the experimental data, Nelson et al. [1993] assumed that the transverse normal stress is equal the vertical normal stress, so k = u2 /2 + v2 . The agreement between simulated and observed turbulent kinetic energy shown in Figure A1 is good for both turbulence models especially in the near wall region, while k is underestimated far away from the bed in the Bennett and Best measurements. Venditti and Bennett [2000] measured all three components of the normal stresses, and they observed that close to the reattachment point the transverse component is significant larger than the vertical normal stress. Farther downstream the magnitude of the two components becomes more or less the same. The k-e model cannot predict this, as it does not produce the individual components in normal stresses. Figure A2 compares the RSM predictions of the normal stresses with measurements. The transverse normal stress prediction is less satisfactory in the separation zone

F04013

Figure A2. Comparison between the RSM model predictions and the Venditti and Bennett [2000] measurements of the three components of the normal stresses normalized by the mean flow velocity. The points in the circle belong to the upstream measuring station as indicated by a thin line.

and upstream part of the wake, while it is satisfactory in the downstream part of the wake and in the boundary layer on the top of the dune, which is the relevant part for the present study. The large transverse normal stresses might have some impact (increase) on the transport capacity of sediment in the trough (see also section 8), and will probably require at least a large eddy simulation flow description.

Notation a c cb cf cf0 cm ce1, ce2 d D g h H Hs k kN L Lm Figure A1. Predicted and measured time-averaged turbulent kinetic energy. Symbols are the same as in Figure 3. (a) Measurements by Nelson et al. [1993]. (b) Measurements by Bennett and Best [1995]. 11 of 13

dune celerity. sediment concentration. sediment concentration at the riverbed. skin friction coefficient. maximum skin friction coefficient. return to isotropy constant. constants appearing in equation (A3). grain diameter. water depth. acceleration of gravity. dune shape. dune height. height of step. turbulent kinetic energy. bed roughness. dune length. location of the maximum bed shear stress downstream a step. n sediment porosity. q sediment transport including pore volume. q0 sediment transport in the dune through.

F04013


qb bed load. qb,top bed load on the dune crest. qD sediment transport deposited behind the dune crest. qs suspended load. qtop sediment transport on the dune crest. s relative sediment density. t time. u horizontal flow velocity. U* friction velocity. ui time-varying flow velocity in the xi -direction. v vertical flow velocity. V depth averaged flow velocity. Vtop depth-averaged flow velocity at the dune top. ws sediment settling velocity. x horizontal coordinate. xi coordinates in the three directions. y vertical coordinate. z Rouse number. a fraction of the suspended load deposited behind the dune crest. dij Kronecker delta. e dissipation of turbulent kinetic energy. k von Karman coefficient. m constant for the inclusion of bed slope in the Shields parameter. nt eddy viscosity. r fluid density. se Schmidt number for the turbulent diffusion of energy dissipation. sk Schmidt number for the turbulent diffusion of turbulent kinetic energy. tb bed shear stress. tk shear stress due to turbulent kinetic energy. tw wall shear stress. q Shields parameter defined in equation (9). q* Shields parameter on a longitudinal slope. qc critical Shields parameter. qtop Shields parameter at the dune crest. qmax iteration value for the Shields parameter. fs angle of repose for sediment. F nondimensional sediment transport. FD nondimensional sediment transport deposited behind the dune crest. Ftop nondimensional sediment transport on the dune crest. [25] Acknowledgment. This research was partially funded by the Danish Research Council (STVF) Research Frame Programs ‘‘Exploitation and Protection of Coastal Zones’’ (EPCOAST), and partly by the EU project HUMOR (contract EVK3-CT-2000-00037).

References Benjamin, T. B. (1959), Shearing flow over a wavy boundary, J. Fluid Mech., 6, 161 – 205. Bennett, S. J., and J. L. Best (1995), Mean flow and turbulence structure of fixed, two-dimensional bed forms: Implications for sediment transport and bed form stability, Sedimentology, 42, 491 – 513. Bradshaw, P., and F. Y. F. Wong (1972), The reattachment and relation of a turbulent shear layer, J. Fluid Mech., 52, 113 – 135. Deigaard, R., and J. Fredsøe (1987), Offshore sand waves, in Proceedings of the 20th Coastal Engineering Conference, edited by B. L. Edge, pp. 1047 – 1061, Am. Soc. of Civ. Eng., Reston, Va. Engelund, F. (1970), Instability of erodible beds, J. Fluid Mech., 42, 225 – 244.

F04013

Engelund, F. (1971), The solitary sand wave, ISVA Prog. Rep. 24, pp. 51 – 54, Tech. Univ. of Denmark, Lyngby, Denmark. Engelund, F., and J. Fredsøe (1976), A sediment transport model for straight alluvial channels, Nord. Hydrol., 7, 293 – 406. Engelund, F., and J. Fredsøe (1982), Sediment ripples and dunes, Annu. Rev. Fluid Mech., 14, 13 – 37. Engelund, F., and E. Hansen (1973), A Monograph on Sediment Transport in Alluvial Streams, Tech. Press, Copenhagen. ¨ ber die Wechselwirkung zwischen Wasser und Exner, F. (1925), U Geschiebe in Flu¨ssen, paper presented at Section IIA, Vienna Acad. of Sci. Fredsøe, J. (1974), On the development of dunes in erodible channels, J. Fluid Mech., 64, 1 – 16. Fredsøe, J. (1979), Unsteady flow in straight alluvial streams: modification of individual dunes, J. Fluid Mech., 91, 497 – 512. Fredsøe, J. (1982), Shape and dimensions of stationary dunes in rivers, J. Hydraul. Eng., 108(8), 932 – 947. Fredsøe, J. (1996), The stability of a sandy riverbed, in Issues and Directions in Hydraulics, T. Nakato and R. Ettema, pp. 99 – 113, A. A. Balkema, Brookfield, Vt. Fredsøe, J., and R. Deigaard (1992), Mechanics of Coastal Sediment Transport, Adv. Ser. Ocean Eng., vol. 3, World Sci., Hackensack, N. J. Guy, H. P., D. B. Simons, and E. V. Richardson (1966), Summary of alluvial channel data from flume experiments, 1956 – 61, U. S. Geol. Surv. Prof. Pap. 462-1, 96 pp. Jensen, J. H., E. O. Madsen, and J. Fredsøe (1999), Oblique flow over dredged channels: II. Sediment transport and morphology, J. Hydraul. Eng., 125(11), 1190 – 1198. Julien, P. Y., and G. J. Klaassen (1995), Sand-dune geometry of large rivers during floods, J. Hydraul. Eng., 121(9), 657 – 663. Kennedy, J. F. (1971), The formation of sediment ripples, dunes and antidunes, Annu. Rev. Fluid Mech., 1, 147 – 168. Klaassen, G. J., H. J. M. Ogink, and L. C. van Rijn (1986), DHL-research on bed forms, resistance to flow and sediment transport, Commun. 362, 25 pp., Delft Hydraul., Delft, Netherlands. Launder, B. E., and D. B. Spalding (1972), Lectures in Mathematical Models of Turbulence, Elsevier, New York. Launder, B. E., G. J. Reece, and W. Rodi (1975), Progress in the development of a Reynolds-stress turbulence closure, J. Fluid Mech., 68, 537 – 566. McLean, S. R. (1990), The stability of ripples and dunes, Earth Sci. Rev., 29, 131 – 144. McLean, S. R., and J. D. Smith (1986), A model for flow over twodimensional bed forms, J. Hydraul. Eng., 112(4), 300 – 317. Mendoza, C., and H. W. Shen (1990), Investigation of turbulent flow over dunes, J. Hydraul. Eng., 116(4), 459 – 477. Meyer-Peter, E., and R. Mu¨ller (1948), Formulas for bed load transport, paper presented at Second Congress, Int. Assoc. for Hydraul. Res., Rotterdam, Netherlands. Mueller, T. J., and J. M. Robertson (1963), A study of the mean motion and turbulence downstream of a roughness element, Mod. Dev. Theor. Appl. Mech., 1, 326 – 340. Nelson, J., and J. D. Smith (1989), Mechanics of flow over ripples and dunes, J. Geophys. Res., 94, 8146 – 8162. Nelson, J. M., S. R. McLean, and S. R. Wolfe (1993), Mean flow and turbulence over two-dimensional bed forms, Water Resour. Res., 29(12), 3935 – 3953. Nelson, J. M., R. L. Shreve, S. R. McLean, and T. G. Drake (1995), Role of near-bed turbulence in bed load transport and bed form mechanics, Water Resour. Res., 31(8), 2071 – 2085. Reynolds, A. J. (1976), A decade’s investigation of the stability of erodible stream beds, Nord. Hydrol., 7, 161 – 180. Richards, K. J. (1980), The formation of ripples and dunes on an erodible bed, J. Fluid Mech., 99, 597 – 618. Richards, K. J., and P. A. Taylor (1981), A numerical model of flow over sand waves in water of finite depth, Geophys. J. R. Astron. Soc., 65, 103 – 128. Rodi, W. (1993), Turbulence models and their application in hydraulics, report, Int. Assoc. for Hydraul. Res., Delft, Netherlands. Smith, J. D. (1970), Stability of a sand bed subjected to shear flow of low Froude numbers, J. Geophys. Res., 75, 5928 – 5940. Sumer, B. M., L. H. C. Chua, N. S. Cheng, and J. Fredsøe (2003), Influence of turbulence on bed load sediment transport, J. Hydraul. Eng., 129(8), 585 – 596. Termes, A. P. P. (1984), Water-movement over a horizontal bed and solitary sand dune. Rep. R/1984/H/8, Dep. of Civ. Eng., Delft Univ. of Technol., Delft, Netherlands. U. S. Army Corps of Engineers (1969), Missouri River channel studies, Mo. River Div. Sed. Ser. 13B, 13 pages, Washington, D. C.

12 of 13

F04013


Van der Knaap, F. C. M., M. C. L. M. Mierlo, and M. J. Officier (1991), Measurements and computations of the turbulent flow field above fixed bed forms, in Euromech 262—Sand Transport in Rivers, Estuaries, and the Sea, edited by R. L. Soulsby and R. Bettess, pp. 179 – 184, A. A. Balkema, Brookfield, Vt. Van Rijn, L. C. (1984), Sediment transport: Part III. Alluvial roughness, J. Hydraul. Eng., 110(12), 1733 – 1754. Venditti, J. G., and S. J. Bennett (2000), Spectral analysis of turbulent flow and suspended sediment transport over fixed dunes, J. Geophys. Res., 105, 22,035 – 22,047. Wilcox, D. C. (1993), Turbulence modeling for CDF, report, DCW Industries, Inc., La Canãda, Calif.

F04013

Yoon, J. Y., and V. C. Patel (1996), Numerical model of turbulent flow over sand dune, J. Hydraul. Eng., 122(1), 10 – 17. Ziljema, M., A. Segal, and P. Wesseling (1995), Finite volume computation of incompressible turbulent flows in general coordinates on staggered grids, Int. J. Numer. Methods Fluids, 20, 621 – 640.

J. Fredsøe, Technical University of Denmark, MEK, Building 403, DK2800 Lungby, Denmark. ([email protected]) S. Tjerry, DHI-Water and Environment, 319 SW Washington, Suite 614, Portland, OR 97204, USA. ([email protected])

13 of 13