Applied Numerical Mathematics Pages 1749–1766.pdf

Applied Numerical Mathematics 62 2012 1749 1766 Contents lists available at SciVerse ScienceDirect Applied Numerical Mathematics Numerical modelling of sediment transport in the Nador lagoon Morocco Fayssal Benkhaldoun a Salah Daoudib Imad Elmahib Mohammed Seaidc aLAGA Universit Paris 13 99 Av J B Clement 93430 Villetaneuse France bENSAO EMSN COSTE Universit Mohammed 1 B P 669 60000 Oujda Morocco cSchool of Engineering and Computing Sciences University of Durham South Road Durham DH1 3LE UK a r t i c l ei n f oa b s t r a c t Article history Available online 28 June 2012 Keywords Nador lagoon Shallow water equations Sediment transport Morphodynamics Finite volume method Unstructured mesh We present a numerical method for solving the sediment transport in the Nador lagoon The lagoon is located on the Moroccan eastern coast and exchanges water fl ow with the Mediterranean sea The governing equations consist of the well established shallow water system including bathymetric forces friction terms coriolis and eddy diffusion stresses To model sediment transport we consider an Exner equation for morphological evolution and an advection diffusion problem for the transport of suspended sediments As a numerical solver we apply an adaptive fi nite volume method using a centred type discretization for the source terms The proposed method can handle complex topography using unstructured grids and satisfi es the conservation property Several numerical results are presented to demonstrate the high resolution of the proposed method and to confi rm its capability to provide accurate and effi cient simulations for sediment transport in the Nador lagoon 2012 IMACS Published by Elsevier B V All rights reserved 1 Introduction The Nador lagoon is located on the Moroccan eastern coast It is a restricted lagoon of 115 km2 25 km by 7 5 km and with a depth not exceeding 8 m see Fig 1 Recently the Nador lagoon has been the subject of many investigations on water quality currents fl ora fauna fi shing and aquaculture see for instance 9 22 Most of these studies deal with the environmental aspects of the lagoon such as biological and socioeconomic impacts However to the best of our knowledge there are no research studies on the numerical modelling of sediment transport in the Nador lagoon Needless to mention that numerical studies are essential since they can quantify the interaction between sediment transport and water fl ow and thereafter can help to understand the evolution of the lagoon morphodynamics Consequently this may provide numerical tools to study the physical environment of the lagoon and to assess the development strategy reducing the fl ood and pollution risks in the lagoon Certainly numerical modelling of sediment transport in the Nador lagoon would be less costly than experimental study on the lagoon fi eld In the current work the governing equations consist of the well established shallow water system including bathymetric forces Coriolis effects friction terms and eddy diffusion stresses To model sediment transport we consider the well known Exner equation for morphological evolution and an advection diffusion problem for the suspended sediment accounting for erosion and deposition effects The coupled hydrodynamical and morphological model forms a hyperbolic system of conservation laws with source terms which are not easy to solve numerically For instance numerical simulation of mor phodynamical changes of the bed in hydraulic systems involve different physical mechanisms propagating within the system Corresponding author E mail address ielmahi ensa ump ma I Elmahi 0168 9274 36 00 2012 IMACS Published by Elsevier B V All rights reserved http dx doi org 10 1016 j apnum 2012 05 010 1750F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 1766 Fig 1 Location and schematic description of the Nador lagoon according to their time response i e the problem is a multi scale system which requires a robust solver to accurately re solve both hydrodynamical and morphodynamical time scales In addition most morphodynamical fl ows frequently involve important features such as internal bores and moving shoreline fronts which present a signifi cant challenge to the accuracy and stability of numerical models The objective of this work is to devise a stable reliable and accurate numerical method able to approximate solutions to the coupled hydrodynamical and morphological system in the Nador lagoon It should be stressed that it is diffi cult to validate numerical results for sediment transport against measurements from fi eld experiments The main diffi culties lie essentially on the empirical formulae needed to close the sediment transport model and also on the calibration of the parameters involved in the bed load modelling we refer the reader to 11 for more discussions An adaptive fi nite volume Non Homogeneous Riemann Solver SRNH has recently been proposed in 6 7 for solving sed iment transport by shallow water equations without accounting for suspended sediment Here the acronym SRNH stands for Solveur de Riemann Non Homog ne The SRNH method belongs to the class of methods that employ only physical fl uxes and averaged states in their formulations To control the local diffusion in the scheme and also to preserve monotonicity a parameter is introduced based on the sign matrix of the fl ux Jacobian The main features of such a SRNH scheme are i the capability to satisfy the conservation property resulting in numerical solutions free from spurious oscillations in signifi cant morphodynamic situations ii the implementation on unstructured meshes allowing for local mesh refi nement during the simulation process and iii the achievement of strong stability for simulations of slowly varying bed load as well as rapidly varying fl ows containing also shocks or discontinuities It should also be stressed that the SRNH method has been success fully applied to solve pollutant dispersion by shallow water fl ows in 3 Neither eddy diffusivity nor wind shear stresses at the water surface have been taken into account in the study presented in 3 Nevertheless our SRNH method has been found to perform effectively on realistic irregular bottoms handle complex geometry to solve different hydraulic regimes and to preserve conservation properties for both water free surface and pollutant concentration Our aim in this paper is to extend the SRNH method for solving the problem of sediment transport in the Nador lagoon We should mention that in the current work we are not interested in solutions with shocks within the Nador lagoon Numerical results presented in this study show that an interesting feature of the SRNH method is to allow multilevel mesh adaptation without deteriorating accuracy of the computed solutions The structure of this paper is as follows In Section 2 we present the mathematical equations for the sediment transport problems considered The formulation of the adaptive fi nite volume method is detailed in Section 3 Section 4 is devoted to numerical results Finally Section 5 contains the conclusions 2 Governing equations for sediment transport problems In this section we describe the physical model used for modelling the sediment transport by shallow water fl ows Here the two dimensional shallow water equations are briefl y recast for the hydrodynamics followed with a short description for the sediment transport equations for the morphodynamics 2 1 Shallow water equations Shallow water equations have been widely used to model free surface fl ows of a fl uid under the infl uence of gravity This class of equations uses the assumption that the vertical scale is much smaller than any typical horizontal scale and can be F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 17661751 Fig 2 Sketch of a domain for two dimensional suspended sediment transport derived from the depth averaged incompressible Navier Stokes equations compare 23 among others For two dimensional fl ow problems these equations are h t hu x hv y 0 1a hu t hu2 1 2gh 2 x huv y gh B x hv bx w wx w Dxx h u v 1b hv t huv x hv2 1 2gh 2 y gh B y hu by w wy w Dyy h u v 1c where u and v are the depth averaged water velocities in x and y direction h the water depth B the bottom topog raphy g the gravitational acceleration wthe water density the Coriolis parameter defi ned by 2 sin with 0 000073 rads 1is the angular velocity of the earth and is the geographic latitude see Fig 2 for an illustration Here bxand by the bed shear stress in the x and y direction respectively defi ned by the depth averaged velocities as bx wCbu u2 v2 by wCbv u2 v2 1d where Cb is the bed friction coeffi cient which may be either constant or estimated as Cb g C2 z where Cz h 1 6 nb is the Chezy constant with nb being the Manning roughness coeffi cient at the bed The surface stress wis usually originated by the shear of the blowing wind and is expressed as a quadratic function of the wind velocity wx wCwwx w2 x w2y by wCwwy w2 x w2y 1e with Cw is the coeffi cient of wind friction and w wx wy Tis the velocity of wind at 10 m above water surface It is usually defi ned by Cw a 0 75 0 067 w2 x w2y 10 3 where ais the air density For the diffusion termsDxx h u v andDyy h u v we have adopted the model derived in 12 16 as Dxx h u v 2 x h 2 u x v y 2 y h u y Dyy h u v 2 x h v x 2 y h u x 2 v y where is the kinematic viscosity Note that we have assumed water fl ow at laminar regime however turbulent effects can also be accounted for in the model 1 by modifying the eddy viscosity coeffi cients In addition other coeffi cients of wind friction can also be applied in 1e 1752F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 1766 2 2 Suspended sediment equations Most of existing formulations of sediment transport are empirical and have been derived based on experimental data and measurements There are many mathematical models for sediment transport by shallow water fl ows In the current work we consider the advection diffusion equation hC t huC x hvC y x h C x y h C y E D hQ 2a where C is the depth averaged concentration of suspended sediment Q is the depth averaged sediment source and is the diffusion coeffi cient of the sediment In Eq 2a E and D represent the erosion and deposition terms in upward and downward directions respectively In the current study the erosion deposition term is calculated as E D s C C 2b where sis the settling velocity of sediment particles is the recovery coeffi cient of the suspended sediment and C is the sediment concentration close to the bed defi ned by 20 C 0 015d50T bd0 3 2c with d50is the mean diameter of the sediment and b is a reference level fi xed in our computations to 0 05h compare 18 The excess bed shear stressT is defi ned as T u2 u 2 c u2 c 2d where the shear velocity u and the critical bed shear velocity ucfor the sediment are given by u2 2 bx 2 by w uc a d 50 10 b with a and b are constants depending on the mean diameter of the sediment particles In 2c the particle size diameter d is defi ned as d d50 g s w w 2 1 3 2e where sis the sediment density Note that the parameters s a and b are user defi ned constants in the sediment transport model 2 3 Bed load equations In 1a the function B corresponds to the sediment layer characterizing the bed level For fi xed bottom topography i e B B x Eq 1a reduces to the standard shallow water equations In the current work we assume that a sediment transport takes place such that the bed level depends on the time variable as well This requires an additional equation for its evolution Here to update the bed load we use the Exner equation given by 1 p B t Qbx x Qby y x B x y B y 3a where p is the sediment porosity assumed to be constant and is the diffusion coeffi cient In Eq 3a Qbxand Qby represent the bed load sediment transport fl uxes in x and y direction respectively These fl uxes depend on the type of sediment and for simplicity in the presentation we consider the basic sediment transport fl uxes proposed by Grass in 13 Qbx Au u2 v2 m 1 2 Qby Av u2 v2 m 1 2 3b where A is a given experimental constant and 1 m 4 is a chosen parameter both of which are specifi c to the particular sediment transport formula In all simulations presented in this paper the parameter m 3 It should be stressed that the fi nite volume method described in this paper can be applied to other forms of sediment transport fl uxes without major conceptual modifi cations For instance the bed load sediment transport functions proposed in 19 21 can also be handled by the proposed fi nite volume method Notice that in practical situations the diffusion coeffi cients and depend on water depth fl ow velocity bottom roughness wind and vertical turbulence compare 15 for more details For the purpose of the present work the problem of the evaluation of diffusion coeffi cients is not considered F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 17661753 Fig 3 Generic control volume and notation For simplicity in presentation we can also reformulate Eqs 1 2 and 3 in a compact conservative form as W t x F W F W y G W G W Q W S W 4 where W is the vector of conserved variables Q is the source term accounting for the slope variations S is the source term accounting for Coriolis forces friction losses and erosion deposition terms F and G are the advective tensor fl uxes F and G are the diffusion tensor fl uxes compare 3 for a similar formulation for shallow water equations over fi xed beds In 4 W h hu hv hC B Q W 0 gh B x gh B y 0 0 S W 0 hv bx w wx w hu by w wy w E D hQ 0 5 F W hu hu2 1 2gh 2 huv huC Qbx 1 p G W hv huv hv2 1 2gh 2 hvC Qby 1 p F W 0 2 h 2 u x v y 2 h v x h C x 1 p B x G W 0 2 h u y 2 h u x 2 v y h C y 1 p B y Note that we have considered only a single suspended sediment with concentration C transported by the shallow water fl ow however the techniques presented in this paper can straightforwardly be extended to sediment transport with multiple species 3 Finite volume non homogeneous Riemann solver The main advantages of the fi nite volume methods lie on their implementation on unstructured triangular meshes and preserving conservation properties of the equations Hence using the control volume depicted in Fig 3 a fi nite volume discretization of 4 yields 1754F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 1766 Wn 1 i Wn i t Ti j N i ij F Wn d t Ti j N i ij F Wn d t Ti Ti Q Wn dV t Ti Ti S Wn dV 6 where N i is the set of neighboring triangles of the cellTi Wn i is an average value of the solution W in the cellTiat time tn Wi 1 Ti Ti WdV where Ti denotes the area ofTiandSiis the surface surrounding the control volumeTi Here ijis the interface between the two control volumesTiandTj nx ny Tdenotes the unit outward normal to the surfaceSi and F W F W nx G W ny F W F W nx G W ny Note that the time stepping in 6 is only fi rst order accurate A second order accuracy in time can be achieved by consider ing a two step Runge Kutta method Its implementation for solving the semi discrete equations 6 is straightforward The fi nite volume discretization 6 is complete once the gradient fl uxesF Wn diffusion fl uxes F Wn and source terms Q Wn and S Wn are reconstructed 3 1 Discretization of gradient fl uxes Let us consider the hyperbolic part in the system 4 W t F W x G W y Q W 7 where the source term Q accounts only of the bed slopes given in 5 Applied to the system 7 the fi nite volume dis cretization 6 yields t Ti hdV Si hunx hvny d 0 t Ti hudV Si hu2 1 2 gh2 nx huvny d gh Si Bnxd t Ti hv dV Si huvnx hv2 1 2 gh2 ny d gh Si Bnyd t Ti hC dV Si huCnx hvCny d 0 t Ti B dV A Si u u2 v2 nx v u2 v2 ny d 0 where 1 1 p Using the local cell outward normal and tangential the above equations can be reformulated as t Ti hdV Si hu d 0 t Ti hudV Si huu 1 2 gh2nx d gh Si Bnxd t Ti hv dV Si hvu 1 2 gh2ny d gh Si Bnyd F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 17661755 t Ti hC dV Si hu Cd 0 t Ti B dV A Si u2 v2 u d 0 8 where the normal velocity u unx vnyand tangential velocity u uny vnx In order to simplify the system 8 we fi rst sum the second equation multiplied by nxto the third equation multiplied by ny then we subtract the third equation multiplied by nxfrom the second equation multiplied by ny The result of these operations is t Ti hdV Si hu d 0 t Ti hu dV Si hu u 1 2 gh2 d gh Si B d t Ti hu dV Si hu u d 0 t Ti hC dV Si hu Cd 0 t Ti B dV A Si u2 v2 u d 0 9 which can be reformulated in differential form as h t hu 0 hu t hu2 1 2 gh2 gh B hu t hu u 0 hC t hu C 0 B t A u u2 u 2 0 10 Hence the system 10 can also be rewritten in a vector form as W t A W W 0 11 where W h hu hu hC B A W 01000 gh u2 2u 00gh u u u u 00 u CC0u 0 A 3u u2 u2 h A 3u2 u2 h A 2u u h 00 Note that one of the advantages in considering the projected system 11 is that no discretization of source terms is required Thus in the predictor stage we use the projected system 11 to compute the averaged states as 1756F Benkhaldoun et al Applied Numerical Mathematics 62 2012 1749 1766 Wn ij 1 2 Wn i Wn j 1 2 sgn A W Wn j W n i 12 where the sign matrix is defi ned as sgn A W R W sgn W R 1 W with W is the diagonal matrix of eigenvalues R W is the right eigenvector matrix and W is the Roe s average state given by W hi hj 2 hi hj 2 ui hi uj hj h i hj x vi hi vj hj h i hj y hi hj 2 ui hi uj hj h i hj y vi hi vj hj h i hj x hi hj 2 Ci hi Cj hj h i hj Bi Bj 2 13 The matrices W andR W can be explicitly expressed using the associated eigenvalues of A W For convenience of the reader these matrices are formulated in Appendix A Once the projected states are calculated in the predictor stage 12 the states Wn ij are recovered by using the trans formations v u u and u u u Thus applied to the system 7 the proposed SRNH scheme consists of a predictor stage and a corrector stage and can be formulated as Wn ij 1 2 Wn i Wn j 1 2 sgn A W Wn j W n i Wn 1 i Wn i t Ti j N i F Wn ij ij ij tQni 14 where Qn i is a consistent discretization of the source te