Numerical. simulation of fluid mud layer in estuaries and coastal areas

Vietnam Journal of Mechanics, NCNST of Vietnam T. XX, 1998, No 2 (1 - 10) ·NUMERICAL. SIMULATION OF FLUID MUD LAYER IN ESTUARIES AND COASTAL AREAS DANG HUU CHUNG, PH.D Institute of Mechanics, 224 Doi Can, Hanoi, Vietnam Tel: 0084-4-8326138 Fax: 0084-4-8333039 Email: dhchnngl@imOl.ac.vn ABSTRACT, In this paper the formation and development of fluid mud layer happening in ·estu~ries and coastal areas are studied in detail through. numerical simulation on the basis of the 2D shallow water equations for tidal flow, the advection-diffusi~n equation for cohesive sediment transport and the equations for fluid mud transport. Numerical solution of a special case for a part of the Severn estuary is obtained using the finite difference method as an illustration of the applicability of the model in practice. On the basis of the resul.ts, initial remarks and evaluation are given. 1. Introduction In,estuaries and coastal regions where there is a high concentration of sediment in suspension, a fluid mud layer is often formed during slack water periods by the process of hindered settling. This amount of sediment comes from the sea or from rivers due 'to the process of flowing through many areas in ·a country or around the sides of mountains. Experiments have established the relationship between settling velocity and cohesive sediment concentration (Tsuruya, Murakami and Irie[l]) and it is seen that a peak value of settling velocity occurs at a concentration of approximately 5kgfm3 . Once the near bed sediment concentration exceeds this value, mud settles towards the bed more quickly than it can dewater and a layer of fluid mud forms. The movement of the fluid mud layer can be described by a restricted form of the shallow water equations. Many complicated physical processes occur at the interface between sediment in suspension and fluid mud and between fluid mud and the rigid bed. These are represented in the model in a parameterised way. Fluid mud can also be formed by waves which can fluidise a muddy bed. However, in this paper, the study is restricted to the case of calm conditions, with a very high cohesive sediment concentration in suspension. These conditions are typical in the Severn estuary during a spring tide. The study is focused 1 on the formulation of the mathematiCal model and the numerical simulation of the model on a computer. The functions describing the processes of exchange between suspended sediment, fluid mud and the bed are explained and the results of the application of the model to the Severn estuary are presented. Quantitative measurements offluidmudjn_J;heS_exenuvere !!_Q_1;_available but the behaviour of the model fits well with the description of fluid mud in the Severn given by Kirby and Parker[2]. It is planned to further develop the model to include the effect of waves. 2. The governing equations and boundary conditions As well known, the most popular assumption used in the mathematical models up to now is that the effect of the bed changing in respect to time on hydrody- namics process is ignored. This is because this effect is insignificant in comparison with the other factors when the average sediment concentration is not large enough. Therefore the mathematical model describing this phenomenon is divided into two separate models: The model of tidal current and the model of sediment transport. The tidal model is the 2D horizontal shallow water equations without the force of wind(Muir Wood and Fleming [3]), fJz + fJ(du) + fJ(dv) = 0 fJt fJx fJy (2.1) (2.2) (2.3) in which x, y are the coordinates, the x axis goes along the shore and the y axis perpendicular to the shore, u and v are the tidal flow velocity components along the x and y directions, respectively, z is the water level above the chart datum, g-acceleration due to gravity, C-Chezy coefficient, d-total water depth, n-Coriolis parameter, D-the eddy viscosity coefficient, and t-time. It should be noted that the terms in the right hand sides of the equations (2.2) and (2.3) present the forces due to the surface slope, rotating of the earth, turbulent diffusion and the friction on the bed in the x and y directions, respectively, the scales of which depend on the situations in question. The equation describing sediment transport is the advection-diffusion equa- tion based on the sediment mass conservation with the exchange between sediment in suspension and fluid mud or mud on the bed taken into account (Roberts[4], Odd and Cooper [5], Odd and Rodger [6], and Le Hir and Kalikow [7]) as follows, 2 8(dc) + 8(qxc) + 8(qyc) = dm + ~(dEx Be)+ ~(dE Be) (2.4) at ax 8y dt ax ax By · y 8y . in which c is the mass suspended .sediment concentration, qx, qy components of discharge per unit of width along the x and y directions, respectively, Ex, Ey the diffusion coefficients for sediment along the x and y, Zb the bed level below the chart datum, Pm the mud density, and dd7 the source-sink term. The model for fluid mud layer includes _the equation of fluid mud mass conser- vation and the equations of momentum conservation that are the restricted form of the 2D horizontal shallow water equations (Roberts [4]), a a a dm at (cmdm) + ax (cmdmum) + By (cmdmvm) = dt (2.5) Bum 1 ( ) Pw 8z t:.p 8zm gdm 8t:.p --+ -- ro-r; -flvm+ -g-+g---+ ---- =0 at dmPm X Pm ax Pm ax 2pm ax (2.6) 8vm 1 ( ) Pw 8z t:.p Bzm gdm 8t:.p -- + -- ro-r; +flum+ -g-+g---+ ---- =0 8t dmPm Y Pm 8y Pm 8y 2pm 8y (2.7) in which Cm is the mass concentration of fluid mud, in general a function of time and space, dm is the fluid mud depth, Urn and Vm are the fluid mud velocity components in the X andy directions, respectively, Zm-the elevation of the interface between fluid mud and sediment in suspension, ro- and r;-the shear stress vectors on the bed and on the interface, respectively, Pw-the water density, and Pm-the fluid mud density, the relationship of which with the water density is the following, Pm = Pw + !:.p, !:.p = 0.62cm. (2.8) From equations (2.6)-(2.7) it should be noted that the external forces making fluid mud move in turn are the shear stress on the bed and on the mud-water interface, the Coriolis force, the slope of water surface, the slope of mud-water interface, and gradient of density that is ignored in the study. About the source-sink term presenting the exchange at" the bed or at the mud- water interface in the equations (2.4)-(2.5), the following processes are introduced: * Erosion: dm ( r ) dt =me Te -1 H(r-- r.), H(x) = { 1, x > 0 0, X::':: 0 {2.9) in which m.(KgfN/s) is the erosion rate parameter, r(Njm2 ) is the actual shear stress at the fluid mud-water interface or at the bed-water interface in the ab- sence of fluid mud, r.-the critical bed shear stress for erosion, and H(x)-the usual Heaviside step function. 3 * Settling of mud from suspension: dm ( r) dt = v8 (c)c 1 ~.rd. H(rd- r), v.(c) = \ VmiJ:u eRa, 1lmin c<-- Ro > Vmin c -- - Ro (2.10) vhere Td is the critical bed shear stress for deposition, v,(m/ s)-the settling velocity, 'min-the minimum settling velocity, and R0 (m4/kg/s) are given from experime~ts. rhis process only occurs on the mud-water interface or on the bed in the case vithout fluid mud. * Entrainment: dm O.l~U dt = VeCmH(lO- R;), Ve = (1 + 63R7} 3/ 4 ' R; = ::~~, ~U2 = (u- Um) 2 + (v- Vm} 2 (2.11) where Ve is the entrainment velocity (m/s), and R;-the bulk Richardson number representing the degree of the flow stratification. Therefore, the entrainment only b.appens when the stratific{'tion of the flow is not strong enough on the water-mud lnterface. * Dewatering (2.12) where v0 is the dewatering velocity (m/s). This phenomenon only occurs on the bed when fluid mud layer exists and the bed shear stress is less than the critical value Td. As mentioned above, the entrainment is a process easily causing the numerical instability, so it requires to treat carefully. To close the problem mathematically the initial and boundary conditions for the situation under consideration are required. They are as follows, * The initial conditions: Due to the area of interest is not large enough, _the initial water surface can be horizontal. The condition for concentration of suspended sediment is subjunctively given so that it should be suitable to the case of mud flow. u(x,y,O) = 0, v(x,y,O) = 0, z(x,y,O) = 12.4 c(x,y,O) = 5, dm(x,y,O) = 0, um(x,y,O) = 0, vm(x,y,O) = 0 (2.13) * The boundary conditions: 4 The 4 kinds of boundaries that are necessary to consider here in turn are the river boundary (x = L), the open sea boundary (x = 0}, the offshore boundary (y = 0} and theJand boundary ((x,y) E Ln}: qt.,(x,y,t)[.,=L = fl(t), q1y(x,y,t)i.,=L= o, c(x,y,t)I.,=L = 5, dm(x,y,t)I.,=L = 0, z(x,y,t)l.,=o = h(t), v(x,y,t)l.,=o = o, (2.14} au I = 0, ae I = 0, v ·_n = o, (x, y) E Ln ax :c=O ax x=O in which q1.,, qty are the components of the total water discharge vector, J;(t) ( i = 1, 2) the given functions at the river boundary, v the flow velocity vector, and n the normal vector unit on the land boundary . . 3. Numerical solutions and a test application Owing to the feature of the model, the equations (2.1}-(2.3} together with given boundary and initial conditions that present the tidal flow have been solved numerically firstly. The ADI method was used, with staggered grid for the deriva- tive in respect to space and Leap-frog scheme for time. The corresponding differ- ence equations are given in previous paper (Chung and.Roberts [8]). The finite difference method is also used to solve the equations (2.4}-(2. 7) together with the initial and boundary conditions (2.13}-(2.14} for the sediment transport mod- el. Specially, QUICKEST (see Leonard [9]} is used for the advection-diffusion equation (2.4} and QUICK [9] for the equation (2.5} to get more accuracy. The difference equations corresponding to the equations (2.4}-(2.7} are the following, (de)~+!= (de)~+ !: [uf_1;e;'_1 - uije;'- (dE.,)i-li (!:) ;_1 + (dE.,);;(!:);] + !~ [viJej- vij_ 1 - (dEy);;(!:);+ (dEY)ij-l (!~);J +t.t(dd7);;' [ f>t ( 1 1 J] n+l _ n f>t 1 n+l n 1 + d" To+ T;., umij - Umij + d" T;.,uii + t.tflvmij mijPm mijPm (3.1} t.t gpw ( n ") f.t gt.p (d" d" d d ) - A - zi+lj- Z;j - ~-- mi+lj- mij + u fPij- u epi+li ' uX Pm uX Pm [ t.t ( 1 1 J] n+l _ n t.t 1 n+l n 1 + dn T0 + T;y Vmij.- Vmij + d" T;yVij - t.tflumij mifPm miiPm (3.2) f.t gpw( n n ) f.t gf.p(d" d" d d ) - A - Z;;- Zij+l - ~ -- mij- mij+l + V epij+l- v ep;j , uy Pm uy Pm 5 dn+l - dn f::.t ( n d* n d* ) Cm miy' - Cm mij + D,.x Cm Umi-lj mi-l- Umij mi + ~~c;;.(v~;jd;,; - v::;ii-1 d;,..T--1) + ~::.t(~7) ,J, i = 2, I- 1, j = 2, J- 1 1 which, n n GRAD;= ci+Ii - cij !::.x CU.RVdi = X { !::.1 2 ( d~i+lj - 2d~ij + d~i-1j)' !::..~2 ( d~i+2j - 2d~i+1i + d~ij)' d;,i = i ( d0· + df:i+l) - f::.t CU RVdf, { !::..~2 ( d~ij - 2d~ij+! + d~ij+2)' CURVc~a = !::..~2 (d~ij-1- 2d~ij + d~ij+l), (3.3) . where udep;j and vdep;j are the bed elevations below the chart datum correspond- ing to the positions of u;i al).d V;j, respectively, and fm-the friction factor on the bed for mud flow. The above difference equation systems are solved over a tidal period. The re- sults of computation at two time points t = 38400s and t = 44400s corresponding to before and after high water, respectively (high water at t = 42000s), are dis- played on Figures (1-4) as the characteristics evaluations over a tidal period. They 6 ' ' ' present the behaviour of the mud velocity field, distributions of fluid mud depths and suspended mud concentration. From here it can be seen that at t = 38400s mud concentration in suspension remains at a level higher than the initial level (5kg/m3 ) nearly everywhere due to the effect of high current, so the fluid mud layer only appears in the area near the land boundary. When t ~ 44400s ·current ' ' -sooo ·0 o!!====;2c;;ol;;ooo~==47:o;i;:o'::'o==;===;;:6o;i;:o'::'o ====;;:so;:i:oc;;o====l Land boundary Fig.l. Contour of fluid mud depth and vector field at t = 38400s -1000 <=' "' "" §2000 0 .c "' " "' :3000 " Q. 0 ·4000 ~ .sooo a!====2;;;o;!,;'oo;;====;;40;!,;'oo;;====;::so~oo:;;=====;;a~ooo:;;;===d Land Boundary Fig. 2. Contour of sus. sed. concentration at t = 38400s 7 -1000 2:' "' , §2000 0 -" "' " U) .:3000 " c. 0 -4000 Land boundary Fig. 9. Contour of fluid mud depth and vector field at t = 44400s -1000 2:' "' , §2000 0 -" "' " "' .:::3000 " c. 0 -4000 Land boundary Fig. 4· Contour of sus. sed. concentration at t = 44400s direction changes first in shallow water, because it has less momentum, while the flow velocity in the middle becomes smaller, so the formation of fluid mud appears and settling flux is large mainly in the area far from the shore. It starts to move 8 with very slow velocity a.s seen from the results, tending to accumulate in the deep channel at slack water. For the cas~ under consideration the maximum mud layer depth is about 0.2m. The peak value of fluid mud layer appears when the water level achieves the minimum value that is explained obviously because of enough small flow velocity. In general the results show a sensible behaviour, compared with the description of mud in the Severn given by Kirby and Parker [5] and show to be applicable in practice. Conclusion A numerical simulation for the formation and moving of fluid mud in estu- aries and coa.stal area.s is implemented and applied to the Severn estuary a.s an illustration example. The finite difference method is used, in which ADI method with staggered grid for the derivative in respect to space and Leap-frog scheme for time is used for tidal flow model. Especially, QUICK and QUICKEST schemes with very high accuracy are applied to the equation of fluid mud ma.ss conservation and the equation of advection-diffusion, respectively in the mud transport model. With the grid size of lOOm x lOOm, the computational area consists of 95 x 50 cells, which is not too coarse and is acceptable in simulation and prediction pur- poses. For fluid mud model, although no data mea.surements are available, but the obtained results show a sensible behaviour, compared with the description of mud in the Severn given by Kirby and Parker [2]. Therefore the model shows to be applicable in practice, and that the numerical simulation can be fully carried out on PC. Acknowledgment This publication is completed with financial support from the Council for Natural Sciences of Vietnam. The author would like to thank Dr W. Roberts for his useful remarks. References 1. Hiroichi Tsuruya, Kazuo Murakami and Isao, Mathematical modelling of mud transport in ports with a multi-layered model-Application to Kumamoto Port- Report of the Port and Harbour Research Institute, Vol.29, No.1, 1990 2. R. Kirby and W. R. Parker, Distribution and behaviour of fine sediment in the Severn estuary and Inner Bristol channel, UK., Canadian Journal of Fisheries and Aquatic Sciences, Vol.40, Supplement Number 1, 1983, pp.83-95 3. A. M. Muir Wood and C. A. Fleming, Coa.stal Hydraulics, Second Edition, Macmilan 4. W. Roberts, Development of a mathematical model of fluid mud in the coa.stal 9 zone. Proc. Instn Civ. Engrs Wat., Marit. & Energy, 1993, 101, Sept., 1'73- 181 5. Odd N. V. M. and Cooper A. J., A twe dimensional model of the movement of fluid mud in a high energy turbid estuary. HR Wallingford, 1988, Jan., Report SR 147 6. Odd N. V. M. and Rodger J. G., An analysis of the behaviour of fluid mud in estuaries. HR Wallingford, 1986, Mar., Report SR 84 7. Le Hir P. and Kalikow N., Balance between turbidity maximum and fluid mud in the Loi:re estuary: lessons of a first mathematical modelling. Int. Symp. On the Transport of Suspended Sediment and its Mathematical Modelling, Florence, 1991 8. Dang Huu Chung and Bill Roberts, Mathematical modelling of siltation on intertidal mudflat in the Severn estuary, Proc. of International Conference on Fluid Engineering, Tokyo, Japan, 13~16 July, 1997, pp. 1713-1718 9. B. P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Computer methods in applied mechanics and engineering 19, 1979, pp. 59-98 .0. R. J. S. Whitehouse and H. J. Williamson, Near-bed cohesive Sediment Pro- cesses - The relative importance of tide and wave influences on bed level change at an intertidal cohesive mudflat site. HR Wallingford, 1996, Feb., Report SR 445 Received June 8, 1998 M6 PRONG s6 LCrP BlJN LONG VUNG cuA soNG vA YEN BIEN Trong ba,i bao nay phm:mg phap mo ph6ng si'l dm;rc Str dvng dEl nghien Ctru Slf hlnh thanh Va phat trign crta lap bun !6ng (t vung etta song Va ven bign mqt each chi tie't. Cac mo hlnh toan hc drrqc stt dvng, d6 Ia cac phrrong trlnh song nrr&c nong 2 chi'eu ngang doi v&i dong tril?m, phrrong trlnh khuye'ch tan doi v&i ntng dq bun cat Ia ltrng va h~ phrr=g trlnh dqng lvc mo ta qua trlnh hlnh thanh va phat tri~n crta bun l6ng. L<ri giai so cho m9t tm(mg hqp cv thg & khu VlfC vung etta song Severn da drrqc thu nh%n nher stt dvng phrrang phap sai phan hi'ru han, drrac xem nhrr mot minh hoa cho tfnh kha thi crta m6 hlnh vao thrrc tiL Du-a • 0 • • • • vao nhiiDg Ht qua nh%n drrqc m<;.t vai nh%n xet danh gia ban dau arrqc neu ra. 10 I 1