Numerical simulation of fluid mud layer under current and waves in estuaries and coastal areas

Vietnam Journal of Mechanics, NCNST of Vietnam T. XX, 1998, No 3 (5- 15) NUMERICAL SIMULATION OF FLUID MUD LAYER UNDER CURRENT AND WAVES IN ESTUARIES AND COASTAL AREAS DANG HUU CH.UNG, PH.D Institute of Mechanics, 224 Doi Can, Hanoi, Vietnam Tel: 0084-4-8326138 Fax: 0084-4-8333039 Email: dhchung@imOLac.vn ABSTRACT. The formation and development of fluid mud layer under the both factors of current and wave that often happens in estuaries and coastal areas are studied in detail through numerical simulation on the basis of the 2D shallow water equations for tidal flow, the advection-diffusion equation for cohes.ive sediment transport and the equations for fluid mud transport. At the same time the model for wave propagation is also included. 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. 1. Introduction Cohesive sediment transport plays an important role in a wide rage of design, maintenance and management problems in estuaries and coastal regions. Accumu- lation of sediment in navigation channels and berths often requires very expensive dredging operations. Once a new development appears, it can cause significant changes to the sediment transport patterns, so a good understanding of the like- ly changes is necessary before any engineering project can proceed. In addition, many pollutants are preferentially adsorbed on to the fine cohesive fraction of the sediment and therefore for environmental reasons it is important to be able to predict the movement of the sediment. 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. Once the near bed sediment concentration exceeds the peak 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 5 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. ali' the other hand, fluid mud can also be formed by waveswhich can fluidise amuddybed;especially, in the case of storm the effect of wave can't be ignored. Here the study is focused on the numerical simulation of the fluid mud layer development with a very high cohesive sediment concentration in suspension on a computer by using finite difference method. The model include 2D shallow water equations for tidal flow, wave propagation equation, the advection-diffusion equation for cohesive sediment transport and the equations for fluid mud transport. A numerical solution for a concrete case is obtained as an illustration example. On the basis of the results, initial remarks and evaluation are given. 2. The governing equations and boundary conditions 1. Tidal model 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 i,gnored. 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 [1]), az + a(du) + a(dv) = 0 at . ax ay (2.1) au au au az . (a2u a2u) ,ju2 + v2 at +u ax+ v ay = -g ax+ Ov + D ax2 + ay2 - gu C 2d (2.2) av av av az (a2v' a2v) ,ju2+v2 at + u ox + v ay = -goy - Ou +.D iJx2 + oy2 - gv C2d (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, 0-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) anq (2.3) present the forces 6 d:c:e to the surface slope, rotating of the earth, t"::·t'b•_:~,~:!.:.~ _:~ff·~~i<JL and the ftiction on the bed in the x and y directions, respectively, the scales of which depend on the situations in question. 2. Wave model The governing equation used by the model can be derived from the time- independent form of the Mild-Slope equation: (2.4) in which TJ is the complex water surface elevation, c-the wave celerity, c9-the wave group velocity, w - the angular wave frequency. By representing TJ = Aeis (A is the wave amplitude and S is the wave phase) and neglecting diffraction effects, the ray tracing model is obtained. In the case when wave height curvature in the main propagation direction (y direction) is much less than in the lateral direction ( x direction), the governing equations are the following: aP aQ (2.5) By ax aQ _.!_(_paQ +kak+~~(_!_a2H)) By - Q 8x By 2 By H 8x2 (2.6) :x(H2MP) + :y(H2MQ) = o (2.7) in which P and Q are components of Wave number vector K, H is wave height, k is the wave separation factor. Wave number is represented in the vector form and K = 'ilS = (P,Q) and Cg M= jKj (2.8) 3. Suspended sediment model 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 [2], Odd and Cooper [3], Odd and Rodger [4], and Le Hir and Kalikow [5]) as follows 8(dc) a(qxc) 8(qyc) - dm a ( ac) a ( ac) at + --ai"- + By - dt + 8x dE. ax + By dEy By (2.9) in which c is the mass suspended sediment concentration, q., qy components of discharge per unit of width along the x and y directions, respectively, E., Ey the 7 diffusion coefficients for sediment along the :t and y, Zb the bed level below the chart datum, Pm the mud density, and ~7 the source-sink term. 4. Fluid mud model 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) + ox (cmdmum) + i)y (cmdmvm) = dt (2.10) GUm 1 ( ) Pw OZ. f::.p OZm dm of::.p --+ ro-r; -rlvm+-g-+g---+g----. =0 at dmPm X Pm ox Pm ox 2pm ax (2.11) OVm . 1 ( ) Pw i)z t:.p OZm dm of::.p -- + ro-Ti +rlum+-g- +g--- +g---- =0 at dmPm · Y Pm oy Pm oy 2pm oy (2.12) in which Cm is the mass concentration of fluid mud, in general a function of time and space, dm is the fluid mud depth, Um and Vm the fluid mud velocity com- ponents in the x ~nd y directions, respectively, Zm the elevation of the interface between fluid mud and sediment in suspension, r0 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 + l::.p, f::.p = 0.62Cm. (2.13) From equations (2.11)-(2.12) it should be noted that the external forces mak- ing 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 (29)-(2.10), the following processes are introduced: * Erosion: dm ( r ) - = m - - 1 H(r - r·) dt e 'Te e ' H(x) = { 1, 0, x>O x~O (2.14) in which m.(KgfN/s) is the erosion rate parameter, r(Nfm2 ) is the actual shear stress at the fluid mud-water interface or at the bed-water interface in the absence of fluid mud, re th~ critical bed shear stress for erosion, and H(x) the usual Heaviside step function. 8 dm ( r) - = v8 (c)c 1-- H(rd- r), dt rd . . v8 (c) = ( Vmin• eRa .. ·vmin c<-- Ro Vmin c> -- - Ro (,, l"' -~· >JJ) where Td is the critical bed shear stress for deposition, v8 (m/ s) the settling velocity, Vmin the minimum settling velocity, and R0 (m4 fkgf s) are given from experiments. This process only occurs on the mud-water interface or on the bed in the case without fluid mud. * Entrainment: (H6) 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 happens when the stratification of the flow is not strong enough on the water-mud interface. * Dewatering (2.17) where v0 is the dewatering velocity (mjs). This phenomenon only occurs on thee bed when fluid mud layer exists and the bed shear stress is less than the critica: 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 cm,c:litions fc.r the situation under consideration are required. They are as follows, * The initial conditions: u(x,y,O) = 0, v(x,y,z) = 0, z(x,y,O) = 12.4, c(x, y, 0) = 5, dm(x, y,O) = 0, um(x, y,O) = 0, Vm(x, y,O) = 0. (2.18) * The boundary conditions: 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 the land boundary ((x, y) E Ln): 9 Field of flow velocities(a) ~--- ---.------ -~ -~-T- !:' 4000 ~-~~-~~-~-~-~--~-~~--~~'::~~~~~~~~~~-~-:~:~:~:~:~~;~~~~~~~~~~~~~~~u ~ ... ------,~~'''''----------------------~~~~ ~ .. '~~~~: ~~~~~~~~~~::::::::::::::::::::~; = 3000 .. ,,,,,,,,,,, ... -------------~-~ 0 .• ,,,,,_, ___________________ _ ,.Q - .................... _____________ , ......... , ·~ :::::::::::::::::::::::::;; m 2000 0 • "' 0 1000 .. ···-.---------------~' I~: m/s I 0~------~------~------~------~--~ 0 2000 4000 6000 8000 Land boundary Contour of Water level(c) 5000 1: ~ I, ) .obl ?_ \'\ ~~ \utd ., ~ (' 3000t ~J\0 \,)d 2000 1000 0 0 2000 4000 6000 8000 10000 Field of flow velocities(b) ;~~~~-~-:-:............ ··-· .. ::-:::::::::· j [: _____ _ i 1000~~~~~~~~~: :; ; •.. :::::::::::: .... 3 3000t'" ~:::::::.:=~~~:~::.··::::::::::::::. _g ~~~;~~~~: ::::::·::. ~ ~~~~~:~- .. ::::::::::: .. !11 2000~ \\, ...... ::::::::·----···"·· ~ ! ' \:~~~~~~~====~~~~~~~~~~::: 1000 1:-12.5 m/s II '\., ------········ ! ~ ~:: ~-- ,.--------r-- 0~--~~~~--~--~~ 0 2000 6000 8000 4000 Land boundary Contour of Water :::l\H) )/)) 3000E- ...,_ r "U " 2000 1000 "t;i ~· ~­ ' level( d) ~ "' ~ \ ... ·~ ~ OL-~--~----~------~----~----~ 0 2000 4000 6000 8000 10000 ~ ~ ..c:: .. ·~ -!:- ~ 0 0 ... ... ... II ~ ~ 0::: " <:::- " 0 0 ... 00 "' II ~ Ol ""' 0 s 0 Ol ..... :-5! ~ ... .a .:: .Si ~ .. ~ " "" s 0 '-' .... 0 ~ ~ ;; ~ " tx: """ ,;, ~ q~x(x,y,t)lx=L = ft(t), qty(x,y,t)lx=L = 0, c(x,y,t)lx=L = 5, dm(x,y,t)lx=L = 0, z(x,y,t)lx=O = h(t), v(x,y,t)lx=O = 0, aul · - -o OX x=O- ' Be I - -o OX x=O- ' v.n =0, (x,y) E Ln, (2.19) in which qtx, qty are the components of the total water discharge vector, f; ( 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 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 derivative in respect to space and Leap-frog scheme for time (Chung and Roberts [6]). For the wave model the Predictor- Corrector method is used, in which the discretisation of eqs. (2.5)-(2.7) is carried out by the row for P, Q and H, respectively. It should be noted that an average oprerator is proposed to improve the solution, that is Q;; = :. (r1Qi-1j-1 + r2Qij-1 + rlQi+li-1 + Q;;), r. = 2r1 + r2 + r (3.1) and fully similar to ( 3.1) for P;; and H;;. The finite difference method is also used to solve the equations (2.9)-(2.12) together with the initial and boundary conditions (2.18)-(2.19) for the sediment transport modeL Specially, QUICKEST (see Leonard [7]) is used for the advec- tion-diffusion equation (2.9) and QUICK [7] for the equation (2.10) to get more accuracy. In order to get an overview for the short term prediction the above differ- ence equation systems are solved over a tidal period, some results of which are illustrated in the figures 1-2. For tidal flow, an important feature of the modec is worth to pay attention, that is the representation of the flooding and drying of the intertidal flats and numerical treatment to ensure stability and accuracy because of the very large tidal range in the Severn. At the same time, the results obtained from tidal model have been compared with the field measurements [8] and discussed in more detail (Chung and Roberts [6]). In this work the results 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 displayed as the characteristics evaluations over a tidal period. It shows that in both the 11 :ase the suspended sediment concentrations are quite high in general. This is :ompletely reasonable, because wave factor actually hinders settling of sediment in suspension and remains concentration at higher value. However, this is only the evaluation for ashort time ofa,ti<fa_Jperi()<fwithout~ig change of wave strength. For the case, in which wave strength is considerably variable, such as before and or---------------------~ ~ ·1000 -1000~ -2000. _f\ -2000 -3000 -~~. ' 1 -3000 ., ,( ' 11' i '!1~1, _,--! ', r -~ -4000 .-4000 I -soooo\====;z;;:oo:;;o;===;_,;!oo:;;o;=='::'so;;';o;;;o;==7.ao;l;o;;o=d. -5000 b=~:;==;:;s;==:;:;:;;:==:;;:>::;;=:::::!. 2000 4000 6000 3000 X X Fig. 2a Fig. 2b Field of fluid mud velocities at t=38400s Contours of mud depth at t=38400s -1000,... -2000- -3000 -4000 1- -5000 l.::==;;:;;;;;;===;::;!;;;:==::;::;:;==::;::t:::::::::::d - 2000 4000 6000 8000 X Fig. 2c Fig. 2d Contours of concentration at t=38400s Field of fluid mud velocities at t=44400s 12 -~6:===,~~~===4~rniD~==~,rniD~===,~~~~ X Fig. 2e Contours of depth at t=44400s -1000 -2000 ·3000 -4000 -5000 6:===;;;:;;;===$:;;::::===;;===~;=:::::!. 2000 4000 8000 3000 X Fig. 2j Contours of concentration at t=44400s after a storm, the fluid mud layer by wave will really become a problem for naviga- tional channels and berths. Also from the Fig. 2 it shows that mud concentration in suspension in the area with fluid mud is smaller than the initial value (4.5kgjm3 in the deep area and 3-4kg/m3 in the shallow area with low currents), while it is higher in the narrow band. As soon as the fluid is formed, it starts to move with very slow velocity as illustrated in Fig. 2, causing serious accumulation i:r the deep channel at slack water. For the case under consideration the maximum mud layer depth is about 0.7m that is really bigger in comparison with the case of no wave (Chung [9]). The peak value of fluid mud layer appears when the water level achieves the minimum value that is explained obviously, because at this time flow velocity is enough small and makes sediment in suspension be able to settle downwards near bed and therefore, fi.uid mud layer depth increases rapidly. Th6 interrelation between suspended sediment concentration and fluid mud depth is analysed and shows suitability on their sensitivity by considering the results c" computation at the two positions corresponding to very shallow and deep areas. In both the cases it shows that when suspended sediment concentration goes down to a minimum value, the mud depth achieves the maximum. And as soon as this drying position becomes flooding again it has the behaviour similar with the oth- er position mentioned above and a similar situation happens somewhere else as shown in Fig. 2. 13 Conclusion The formation and development of fluid mud layer in estuaries and coastal areas are numerically simulated and successfully applied to the Severn estuary as an illustration example. Using the finite difference method, numerical solutions for the mathematical models are obtained, 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. For wave model the predictor-corrector method shows an effect. Especially, QUICK and QUICKEST schemes with very high accuracy are applied to the equation of fluid mud mass 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 purposes. Computed results are quite ~'Sensible and show the feasibility of the models. Acknowledgment This publication is completed with financial support from the National Basic Research Program in Natural Sciences. The author would like to thank Dr Bill Roberts for his useful remarks. References 1. Muir Wood A. M. and Fleming C. A. Coastal Hydraulics, Second Edition, ·Macmilan 2. Roberts W. Development of a mathematical model of fluid mud in the coastal zone, Proc. Instn Civ. Engrs Wat., Marit. & Energy, 101, 1993, 173-181. 3. Odd N. V. M. and Cooper A. J., A two dimensional model of the movement of fluid mud in a high energy turbid estuary, HR Wallingford, 1988, Jan., Report SR 147. 4. 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. 5. Le Hir P. and Kalikow N. Balance between turbidity maximum and fluid mud in the Loire estuary: lessons of a first mathematical modelling, Int. Symp. On the Transport of Suspended Sediment and its Mathematical Modelling, Florence, 1991. 6. 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. 14 7. 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 8. Whitehouse R. J. S. and Williamson H. J. 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. 9. Dang Huu Chung. Numerical simulation of fluid mud layer in estuaries and coastal areas, (submitted to J. of Mechanics, 1998) Received June 22, 1998 MO PRONG so LOP BlJN LONG VUNG ctr A SONG VA VEN BIEN DUCJI TAC ElQNG Cl'JA DONG VA SONG sv hlnh thanh va phat trign 16-p bun 16ng (y vung ell-a song va ven bi~ dtr&i Slf tac d{ing cUa d. hai yiJu to dong va song da dtrqc nghien CUou chi ti~t thOng qua phtrang phap mo ph6ng so. Mo hlnh tmin hQC mo t<l. qua trlnh v~t ly phrrc t~p nay bao g'Om h~ phtrang trlnh ntr&c nong 2 chih ngang doi v&i dong tri~u, phtrang trlnh khuyiJch tan xac djnh nang d{i bun cat lo- ltrng va. h~ d.c phtrang trlnh mo ta Slf phat trign ella 16-p bun 16ng. D'Ong thai h~ phtrang trlnh doi v&i khuc X~ song vung ntr&c nong ciing dtrqc str d~ng. LCri giru so cho 1 tmerng hqp c~ thg & vung ell-a song Severn b~ng phtrang phap sai pharr hi'ru h~ da thu nh~n xem nhtr m{it fuinh hQa cho vi~c ap d~ng mo hlnh trong vi~c giai quyiJt cac bai toan trong thvc t~. 15