Coastal Engineering, 15 ( 1991 ) 247-256
247
Elsevier Science Publishers B.V., Amsterdam
A model for sediment re-entrainment and transport in shallow basin flows Karl J. Eidsvik and Torbjorn Utnes Norwegian Hydrotechnical Laboratory, Klaebuveien 153, N- 7034, Trondheim-NTH, Norway (Received December 12, 1989; accepted after revision August 22, 1990)
ABSTRACT Eidsvik, K.J. and Utnes, T., 1991. A model for sediment re-entrainment and transport in shallow basin flows. CoastalEng., 15: 247-256. A finite-element vertically-integrated mean flow model and a boundary layer turbulence/sediment model are combined. The resulting model is "first-order" three-dimensional, enabling predictions of sediment re-entrainment, transport and fallout in shallow basins of arbitrary, smooth shape. The model is well behaved and appears to be realistic.
INTRODUCTION
Sediment re-entrainment is believed to be important for sediment transport. This is, for instance, one of the main reasons why significant parts of the sediments released by European rivers north of Belgium could eventually deposit in the Norwegian Trench (Siindeman and Puls, 1988). This study addresses modelling of sediment re-entrainment and transport in shallow rivers, lakes, oceans, harbours and sediment traps. The flow is supposed to be characterized by vertical shear layers. The geometry of the basin is supposed to be important, while "secondary flows" are considered not to be important. Many aspects of such flows have traditionally been predicted accurately with vertically-integrated models. However, when investigating sediment re-entrainment and transport, the bottom stress, turbulent structure and vertical profiles should be described in detail. A three-dimensional model with turbulence and bottom interaction is therefore needed (Wang, 1989, with references therein). There are several reasons why the three-dimensional aspect of a model should be minimized. Even the most efficient computers can only simulate three-dimensional flows over a grid of the order of 603. This may be sufficient for expected value flow estimations (Rodi, 1987). However, even if a full three-dimensional expected value model is conceivable over a realistic do0378-3839/91/$03.50
© 1991 - - Elsevier Science Publishers B.V.
248
K.J. EIDSVIK AND T. UTNES
main, it is very expensive. Due to its large size, it may also be uncertain. Instead, the model should be tailored towards the difference in horizontal and vertical scales. The tailoring should also be done with the main purpose in mind, which is the estimation of boundary layer structure and sediment reentrainment. With large difference in horizontal and vertical scales, as in the present case, an economic grid of ca. 202 × 60 should be sufficient. This is an order of magnitude smaller than that necessary for a fully three-dimensional model. High peak values for bottom stress, essential for sediment entrainment, may be effectively generated by waves (Hagatun and Eidsvik, 1986; Battjes, 1988 ). Bedload sediment transport is expected to be important when the bottom stress is marginally larger than the critical stress (Graf, 1971 ). Nevertheless, in this study only stress imposed by wind-generated turbulence and suspended sediments are considered. With sediments, the flow is stratified with possible mean velocity maxima. As illustrated by Brors and Eidsvik (1989), these features may require dynamical modelling of Reynolds stress and scalar fluxes. However, the simpler standard (k-e) model without stratification terms, is also expected to be reasonably realistic (Eidsvik and Brors, 1990). For the illustrative purpose of this study the latter is applied. To economize representation of the slow horizontal variations and simultaneously enable nontrivial basin geometry, the model for horizontal flow variations is posed in finite-element formulation by the use of a verticallyintegrated model (Utnes, 1989). The pressure gradient from this model is applied as forcing for the turbulence model. The bottom stress from the turbulence model is applied as forcing in the vertically integrated model. The combined model is "first-order" three-dimensional. It turns out to be well behaved with specific predictions. The predictions appear to be as realistic as the data allow us to decide. MODEL Vertically-integrated model The model is formulated in terms of vertically-averaged, expected value horizontal velocity ai ( i = 1, 2) and surface elevation q as follows (Utnes, 1990): Oai O~i ~ Orl [-l~j~-d-~-i3jJ3Uj=--g + [fi(h+rl) ] -1 (zsi- Zb,) u.~j Ot #Xi Oq+ 0
(1)
(h+,) aj=0
Here the time and horizontal Cartesian coordinates are t, Xl, x2, respectively.
SEDIMENT RE-ENTRAINMENT AND TRANSPORT IN SHALLOW BASIN FLOWS
249
The Coriolis component, gravity acceleration, density, depth and surface elevation are: f3, g, P, h, r/. The wind stress at the surface is zsi. In conventional models of this type the bottom stress, zbi, is set proportional to the verticallyintegrated velocity or its square. Jamart and Ozer ( 1987 ) have discussed deficiencies of such formulations. A simple feature as a bottom Ekman layer, with nonparallel stress and free flow, is sufficient to reveal that more realistic closures are desirable. In the present study Zbi is obtained from the dynamic turbulence model. Boundary conditions at the basin boundary are specified in terms of mean velocity and normal derivative of surface elevation: ui=0; Oq/OXn=Tsn [ p g ( h + q ) l -~ . The time integration is based upon a time-splitting algorithm. The space discretization is performed by a finite-element method with biquadratic elements, and the interpolation is mixed with bilinear elevation and biquadratic velocity. This method is described in detail by Utnes ( 1990 ).
Turbulence- sediment model Vertical integrated models have also been applied to turbulence and sediment characteristics (Pavlovic et al., 1987 ). However, since the vertical sediment variations are expected to be essential, a different approach is chosen here: The slow horizontal variations of the flow means that turbulence and sediment characteristics do locally appear as for a thin shear layer. This implies that a boundary layer turbulence-sediment model, with only horizontal advection terms, should be appropriate. With boundary layer normal coordinate z, and horizontal coordinates (x~, x2) such a model is formulated as:
Ou, ~_u OUi+ 10p + O(vOUi) Ot OXj e'3J3uj=-pOx---~ic°s°~'+(p-pf)p-m gsin°ti Oz\ ~z ]
Oc Oc- o (cw~)+o (~ Oc)
~+uj
Oxj Oz
Oz ~ Oz
(2) Ok
~+
Ok
OuiOu i
Oxj
Oz Oz
u~--=.--
ofvok\
-~+~l,
C~ v Oz Oz
Oz \a, Oz/
The indices and summation convention apply to the horizontal components only. The standard ( k - e ) turbulence model without gravity terms is applied (Eidsvik and Brors, 1989 ). The state variables (ui, c, k, e ) are expected value horizontal velocity components, sediment concentration, turbulent kinetic
250
K.J. E1DSVIK AND T. UTNES
energy and turbulent dissipation. The expected density is p = (1 --c)pf"F COs, with fluid and sediment density pf and Ps, respectively. The sediment fall velocity is ws and the turbulent viscosity coefficient is u= C~, k2/e. Since the turbulence is supposed to be locally controlled, with the most significant production near the surface and bottom, the coordinates (x~, x2, z) may also be considered to be weakly curvilinear in a boundary fitted frame. The particular pressure gradient-gravity forcing term in the momentum equation originates from transformation of the pressure gradient to a coordinate system aligned with the bottom, with local slope angle ai. Along the vertical axis the hydrostatic and Boussinesq approximation is applied: Op/Oz~ -gpf. Near the smoothly varying bottom this should enable a firstorder approximation to simulate turbidity forcing. Near the surface this term vanishes because p--.pf there. In the shallow basin, the horizontal pressure gradient has small vertical variations. Therefore Op/Oxi in Eq. 2 is approximately equal to Op/Oxi=pg Orl/Ox~ in the vertically-integrated model. Experimental coefficients are (C,, C~, C2, ac, Crk,or,) = (0.09, 1.44, 1.92, 1.0, 1.0, 1.3 ). These particular values are chosen because they are the most commonly used. Obviously missing from Eq. 2 is the continuity equation and the vertical convective terms. These aspects impose horizontal coupling and vertical velocity field. For the purpose of studying the simplest possible model, the vertical velocity and horizontal diffusive terms are assumed to be negligible. Continuity is imposed only in the vertically averaged sense: At each time-step the profile, ui(z), is adjusted additively, so that the correct integral value is obtained:
u,(z),-
(3)
Boundary conditions at the surface are specified as an equilibrium boundary layer. With large uncertainty on the structure of the oceanic boundary layer towards the atmosphere (Gargett, 1989), a simple alternative with a surface roughness OfZo = 10-3 m is chosen. At the bottom, Hagatun and Eidsvik's (1986) boundary conditions are applied, with U(Zo) = 0 at Zo= 10 -3 m. The vertical grid is chosen to be consistent with the weakly curvi-linear system mentioned. Near the surface and bottom the vertical grid is refined as geometrical series. With a constant number of vertical grid points, the grid is stretched in deep water and refined in shallow water. INTEGRATIONS
It is difficult to be certain that a model with so many-dimensional stateand parameter-spaces and so few data really has a unique and realistic solution. It is even a point to ensure that serious errors are avoided. However, each of the two models have separately been indicated to give well defined and realistic integrations (Eidsvik and Brors, 1989; Utnes, 1990). The prop-
SEDIMENT RE-ENTRAINMENT AND TRANSPORT 1N SHALLOW BASIN FLOWS
251
erties of the combined model are sampled by simulations of flows for which previous experience exists.
Convergence properties, data comparison The question of convergence is investigated by estimating the approach to a stationary flow applying a different vertical grid resolution. The flow chosen for this study is a wind-forced flow in a rectangular basin with dimensions 200 m alongwind and 100 m crosswind. The constant wind stress is chosen as 9.2X 10 -2 ( N / m 2) equivalent to a surface friction velocity u , = 9 . 6 X 10 - 3 m/s. The time-step is chosen as 30 s. Figure 1 illustrates the time development of the basin center velocity profiles (towards stationarity?), with and without correcting for vertically-integrated continuity. Without correction there is at best a large time-scale error. This is sufficient reason for imposing vertically-integrated continuity. The corrected flow model converges within a few time-steps (in shallow water without Coriolis terms). The vertical grid dependence is illustrated in Fig. 2. The number of vertical grids is varied through 30, 50 and 100. As observed, the mean velocity profile does not reveal significant grid dependence. However, the turbulent viscosity profile shows slight grid variations in the lower half of the basin. Whether the integrations converge towards the correct result, is judged by comparisons with data from Baines and Knapp ( 1965 ) and predictions by Svensson ( 1978 ) and Sheng (1984). As seen from Fig. 2 there is a remarkable similarity be-
II //
/ /
I
iI iI
-s
~
~
~'o u/o,
~s
Fig. 1. Time development towards a stationary wind forced flow in the middle of a rectangular basin with alongwind and crosswind dimensions of 200 and 100 m, resp. The number of timesteps is represented by 100, 200 and 400. Time-step: A t = 30 s; - - - Noncorrected turbulence model; Vertically integrated continuity imposed, with convergence before 200 time-steps.
252
K.J. EIDSV1K AND T. UTNES
1.0
z/h .8 ~ \ \ \ 150 ~'-~h
!
l I
-5
i
II
//
i
0
15
Fig. 2. Grid variations of stationary wind forced flow in the middle of a rectangular basin. Only the turbulent viscosity coefficient is significantly affected. The number of vertical grids is represented by 30, 50 and 100. Open circles represent data by Baines and Knapp ( 1965 ). tween the data and the predictions, without parameter adjustments o f any kind. Close to the surface and bottom, where measurements are difficult, the present model predicts essentially the same m e a n velocity profiles as Svensson (1978) and Sheng (1984). The m i n i m u m u/u, is approximately - 2 . 5 at z/h = 0.1 in all predictions. Near the surface the m e a n flow gradients are so large and sensitive to small variations in surface mixing, that it is difficult to compare models. It appears that Svensson's and Sheng's surface values for u/ u, are approximately 20 and 10, respectively. The present Figs. 1 and 2 suggest a value o f the same magnitude. Also the turbulent viscosity profiles turn out to be predicted similar, with m a x i m u m u/u,h = 8 X 10 -2 at z/h = 0.55. The predictions illustrated support a general experience that the model is numerically well behaved with rapid convergence and limited vertical grid dependence provided that the n u m b e r o f grid is of the order o f 100. The integrations appear to be reasonable as far as our data and intuition reach.
Flow in a circular basin To exploit the potential for simulating the flow in realistic basins with arbitrary smooth shape, a circular basin with parabolic depth distribution is chosen as an example. The horizontal basin size and grid are illustrated in Fig. 3. Time convergence is found to be rapid also in this flow in that the solution changes insignificantly from one to two days integration time. Figure 3 shows the predicted flow without convective terms. Integrations with convective
SEDIMENT RE-ENTRAINMENT
8
AND TRANSPORT
253
IN SHALLOW BASIN FLOWS
(a)
~_ (b)
E
0.5 m / s
Grid
8 8
//
//
\
\
rl
\,
\
I
--
8
//
\ /
,-%\\
o
//
l
o I
-5000
-~500
2500
2500
5000
m
7500
8
(c)
-5000
-7500
75O0 (m)
(d) 0.5 m / s I
ui ("q- a z )
I
\ '\ \.\ ~, ~,
10
\
:/'/1
i
~//
'\
' \ '\)ti\
I
/ /
: :/IL 8
'\ % . \ C:)
\
!//
o
-~500
-5000
-2500
0
2~00
5000 (
m
~. (e) o
k(Az)
:////
~\%,u.
) 7500
8
(fi
8
I
":,
-2~
0.5m/s
(m)
I
ui [Az)
8
<%\
8 ; ~,,%
-5~
-7500
,,,,,,
/
/
T
L"Ix"'i'
8 8 8
, \
\ \
\
\
I
I
/
o
(m)
-7500
i
Fig. 3. Stationary wind forced flow in a circular basin on rotating earth. Parabolic depth distribution with ( m i n i m u m , m a x i m u m ) : (3, 10) m. Time-step At= 5 rain; Integration time: 48 h; Southerly wind stress: 0.1 N / m --2. Isocurves for q: ( X 10 -3 m ) , for k: ( X 10 -5 m 2 s-2); ~Jz: Grid difference.
254 8
o!P
K.J. EIDSVIKAND T. UTNES (a) i
(b) 0.5 m / s
k
u i (&z)
f
~-t \\
v 4
o
.
g,
x
/
t
/
o ~7500
(&Z)
o '
r
5000
2500
0
[
I
2500
5000
I
(rn)
75O0
7500
~ 2500
0
2500
5000
(m)7500
~7 (c) c (lOAz)
c (L~z)
g
< /
w
oo
8
o
o
o r
-7500
-5000
i
-2500
0
25o0
5 6 o 0 ( m ) 7 5io o
-750o
-5000
-2500
0
2 oo
5000 (m) 7 oo
Fig. 4. Wind forced flow in a circular basin. Parameters as in Fig. 3 except: Convective terms included. Depth ( m i n i m u m , m a x i m u m ) : (1.5, 5) m; Integration time: 24 h; Southerly wind stress: 0.2 N/m -2. lsocurves for k: ( X 10 -5 s-2), for c: ( X 10-2); Az: Grid difference•
terms included give similar predictions (Fig. 4 ). Predictions without the Coriolis term (not shown) are also similar, but without the asymmetry evident in Figs. 3 and 4. Giorgini and Gray (1980) and Liable (1984) have previously simulated the flow in such a basin, using algebraic turbulence models. Their predictions are similar to ours, suggesting that no serious errors are made and that wellbehaved predictions are obtained. The vertically-integrated velocity field, representative for the mass transport, shows that the alongwind mass transport occurs near the sides of the basin. In the central portion of the basin, there is a mass transport in the opposite direction of the wind stress. The surface elevation variation is small,
SEDIMENT RE-ENTRAINMENT AND TRANSPORT IN SHALLOW BASIN FLOWS
255
with the strongest surface currents near the sides. The central basin surface current and elevation approach geostrophic balance. The bottom current and kinetic energy are also reasonable. Since the height of the lowest grid point is larger in the deeper basin center than near the shallower basin sides, the magnitude of the bottom flow is exaggerated by Fig. 3 near the basin centre. The turbulent kinetic energy is insignificantly affected by this. Since the absolute value of the bottom stress is: Izl =PC 1/2 k, with direction given by the mean velocity close to the bottom, Fig. 3 also describes the bottom stress vector. In the central portion of the basin there is a systematic counter-wind flow with weak bottom stress. Near the basin sides there is a relatively strong alongwind flow with larger bottom stress. The sediment particles are assumed to be cohesionless with diameter d----5X 10 -4 m. Based upon a Shields number criterium for entrainment, S = Izl [ (p~-pf)gd] - ~>~0.1, the bottom stress in Fig. 3 is not large enough to entrain these sediments. Figure 4 illustrates the predicted bottom flow in a shallower basin with high enough wind stress for entrainment. The Shields number is still small, and the sediment concentration field turns out to change significantly due to small variations of the integration method. The details of the concentration asymetry shown in Fig. 4 may therefore be unrealistic. However, the main feature of m a x i m u m concentrations near the shallow basin sides appear to be realistic. With bottom stress and sediment type given, a standard bedload transport model (compare Graf, 1971) can easily be implemented. CONCLUDING REMARKS
A finite-element vertically-integrated mean flow model and a boundary layer turbulence/sediment model are combined. The result is a first-order threedimensional model enabling predictions of sediment re-entrainment, transport and fallout in basins of arbitrary smooth shape. The model is well behaved. It predicts cheaply and appears to be realistic. ACKNOWLEDGEMENT
The study was funded by The Fund of Licence Fees and Norwegian Hydrotechnical Laboratory.
REFERENCES Baines, W.D. and Knapp, D.J., 1965. Wind-driven water currents. J. Hydrol. Div. ASCE., 91: 205-221. Battjes, J.A., 1988. Zurf-zone dynamics. Ann. Rev. Fluid Mech., 21: 275-293.
256
K.J. EIDSVIK AND T. UTNES
Brors, B. and Eidsvik, K.J., 1989. Selfaccelerated turbidity current prediction based upon dynamic turbulence models. Preprints: ASCE Int. Symp. Sediment Transport Modelling. New Orleans, L.A. Davies, A.A. and Jones, J.E., 1988. Modelling turbulence in shallow sea regions. In: J.C.J. Nihoul and M.M. Jamart (Editors), Small-Scale Turbulence and Mixing in the ocean. Elsevier, Amsterdam, 542 pp. Eidsvik, K.J. and Brors, B., 1989. Selfaccelerated turbidity current prediction based upon ( k - ~ ) turbulence. Cont. Shelf Res., 9:617-627. Gargett, A., 1989. Ocean turbulence. Ann. Rev, Fluid Mech., 21: 419-451. Giorgini, A. and Gray, D.D., 1980. Wind-driven circulation in small shallow lakes: Laboratory and computational experiments. Purdue Univ. Water Resour. Res. Center, West Lafayette, Ind., 103 pp. Graf, H.W., 1971. Hydraulics of Sediment Transport. McGraw-Hill, New York, N.Y., 513 pp. Hagatun, K. and Eidsvik, K.J., 1986. Oscillating turbulent boundary layer with suspended sediments. J. Geophys. Res., 91: 13045-13055. Jamart, M.M. and Ozer, J., 1987. Comparison of 2-D and 3-D models of the steady wind-driven circulation in shallow waters. Coastal Eng., 11 (5/6): 393-413. Liable, J.P., 1984. A finite element/finite difference wave model for depth varying nearly horizontal flow. Adv. Water Resour., 7: 2-14. Pavlovic, R.N., Varga, S. and Misic, B., 1987. Two-dimensional depth-averaged model for the calculation of sediment transport and river bed deformation. In: C.J. Chen and F.M. Holly (Editors), Turbulence Measurements and Flow Modelling. Springer, Berlin. Rodi, W., 1987. Examples of calculation methods for flow and mixing in stratified fluids. J. Geophys. Res., 92: 5305-5328. Sheng, P.Y., 1984. Turbulent transport model of coastal processes. In: Proc. ASCE, 19th Coastal Eng. Conf., pp. 2380-2396. Siindemann, J. and Puls, W., 1988. A numerical model of suspended sediment transport in the North Sea. In: 2nd Int. Symp. Wave Research and Coastal Engineering. Svensson, U., 1978. A mathematical model for the seasonal thermocline. Ph.D. thesis. Imperial College, London. Utnes, T., 1990. A finite element solution of the shallow water wave equations. Appl. Math. Modelling, 14: 20-29. Wang, S.S.Y., 1989. Preprints: ASCE Int. Symp. Sediment Transport Modeling. New Orleans, L.A., 829 pp.