Table of Contents Author Guidelines Submit a Manuscript
Modelling and Simulation in Engineering
Volume 2014 (2014), Article ID 179504, 8 pages
Research Article

Numerical Survey of Contaminant Transport and Self-Cleansing of Water in Nador Lagoon, Morocco

1LME, Faculté des Sciences, BP 717, 60000 Oujda, Morocco
2School of Engineering, The University of Edinburgh, The King’s Buildings, Mayfield Road, Edinburgh EH9 3JL, UK
3ENSAO, EMCS, Complexe Universitaire, BP 669, 60000 Oujda, Morocco

Received 18 January 2014; Revised 10 May 2014; Accepted 12 May 2014; Published 27 May 2014

Academic Editor: ShengKai Yu

Copyright © 2014 E. M. Chaabelasri et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Numerical simulations are presented of the flow hydrodynamics and hypothetical contaminant dispersion patterns in Nador Lagoon, a shallow lagoon with a barrier island situated on the coast of Morocco. It is found that the natural circulation forced by the tidal flow in the lagoon is greatly affected by the development of an artificial inlet in the barrier island. The case study demonstrates the potential use of modern computational hydraulics as a tool integrated in the decision support system designed to manage a lagoon ecosystem.

1. Introduction

The Nador Lagoon (Figure 1) is located on the eastern coast of Morocco. It is a barrier lagoon of plan area about 115 km2 and horizontal dimensions about 25 km by 7.5 km. The lagoon is shallow with a water depth not exceeding 8 m. Its barrier island is 25 km long with an average width of 300–400 m (reaching 2 km in certain places) and is cut by an artificial inlet. The lagoon is nourished by several tributaries which nowadays serve as receiving waters for sewage discharged from upstream outfalls. There are three main urban centres in the surrounding area (which contribute to the sewage discharge) and two saline workings on the banks of the lagoon. The lagoon also hosts a fish farm that is of considerable importance to the local economy.

Figure 1: Map of the Nador Lagoon with inset view of the bathymetry contours and showing the triangular mesh used by the numerical model.

The lagoon ecosystem is particularly sensitive to water contamination; and its evaluation and analysis are vitally important from both environmental and socioeconomic points of view. Reliable and computationally efficient estimates of the likely deterioration in water quality due to contamination could inform government regulations concerning environmental protection of the lagoon. Moreover, the implementation of counter measures is impossible without an accurate picture of various key processes, including the stage-discharge inflow conditions from the feeder river system, the transport of contaminants within the lagoon itself, and its capacity for self-cleansing. To this end, computational shallow flow simulations supplemented by species advection-diffusion models (representing the transport, dispersion, and fate of contaminants) are useful both for the design of long-term prevention and mitigation measures and for the day to day management of the lagoon. Although no field data are presently available for the Nador Lagoon, several nonvalidated theoretical studies have been undertaken into contaminant and bed-load transport in the lagoon [1, 2].

In this study, a two-dimensional finite volume model is constructed of the Nador lagoon, which couples the shallow water equations and the species transport equation to simulate the transport of passive contaminants by tidal and river circulation within the lagoon. The resulting coupled equations are expressed as a hyperbolic system of conservation laws with source terms [3, 4] and discretized spatially using a finite volume scheme [5] such that mass and momentum are conserved in each computational cell, even in the presence of flow discontinuities. A Godunov-type shock-capturing scheme is used to provide high resolution of the nonhomogenous hyperbolic system of equations on a triangular unstructured computational mesh that is boundary fitted to the Nador Lagoon. Roe’s approximate Riemann solver is used to determine the advective fluxes at cell interfaces [6]. Time integration is by means of a second order Runge-Kutta algorithm.

The main advantages of the finite volume method used herein are (i) its implementation on unstructured meshes that fit a complicated 2D domain; (ii) the simultaneous advection in time of the water flow and the pollutant concentration, solving both problems with the same order of accuracy; (iii) the ability to handle calculations of slowly varying flows or concentrations as well as rapidly varying flows containing shocks and/or fronts; (iv) satisfaction of the exact C-property; and (v) guaranteed positivity of water levels and pollutant concentrations in the transient simulations. All computations are undertaken on a personal computer with Core 2 Duo processor of 2.10 GHz CPU.

This short paper is structured as follows. Section 2 states the governing shallow water and species transport equations. Section 3 briefly describes the numerical solver. The Nador Lagoon simulations are presented in Section 3, and the results are interpreted in terms of the impact of an artificial inlet on the lagoon circulation system and its flushing characteristics. The results provide insight into the likely behavior of pollutants in the Nador Lagoon and indicate hydrodynamic flow structures due to its complicated bed topography and friction. Finally, the water quality renewal capacity of the Nador Lagoon is investigated by considering the self-cleansing processes for two artificial inlet options. Section 4 concludes the paper.

2. Mathematical Model

For shallow flow domains, such as the Nador Lagoon, where the flow is mainly horizontal, the vertical acceleration can be ignored and hydrostatic pressure is assumed. This implies that all waves simulated are long waves, whose amplitude is much smaller than both the depth and the wave length. Moreover, the domain is sufficiently small that the effect of Coriolis acceleration owing to the Earth’s rotation can also be ignored. In such cases, an adequate mathematical description of the flow hydrodynamics is provided by the two-dimensional (2D) shallow water equations. These are coupled with advection-diffusion equation to represent the transport of waterborne contaminants. The resulting equation system is written in tensor notation as [7] where and are indices and the Einstein summation convention is used; that is, repeated indices mean a summation over the space coordinates; is the Cartesian coordinate; is water depth; is time; is the depth-averaged velocity component in the th direction; is the bed elevation above a fixed horizontal datum; is gravitational acceleration; is water density, is depth-averaged concentration, is the dispersion coefficient in the th direction, is the depth-averaged source term; is bed shear stress in the th direction defined by in which is the bed friction coefficient, which may be either constant or estimated from , where is the Chézy constant ( is the Manning coefficient) and is the wind stress components defined by with the coefficient of wind friction defined by where is the air density and is the velocity of wind in the th direction.

Equations in (1) are written in the following conservative form: where and are the vectors of conserved variables and source terms, is the advective flux vector, and is the diffusion flux vector

3. Numerical Formulation

The integral of (2) over a triangular element is written as the sum of each edge contribution where is the vector of conserved variables evaluated at time level , is the number of time steps, is the time step, and is the cell area . To evaluate the state , an approximation is required of the convective and diffusive flux terms at each cell edge. The integral along the edge of a control volume of the normal flux can be written as where is the numerical flux vector and is the edge length of . Herein, the following upwind scheme based on Roe’s approximate Riemann solver is employed to determine on the control volume surfaces. At each cell edge and at time [8], in which where is the flux Jacobian evaluated using Roe’s average state, and are the right and left eigenvector matrices of , and is a diagonal matrix of the absolute values of the eigenvector of . A four-point finite volume interpolation is used to evaluate the diffusion flux through an inner edge , so that where and is the intersection of the orthogonal bisectors of the edge , is the distance between and the edge , and is the length of the edge.

The source terms are balanced by means of a two-dimensional implementation of the upwind scheme proposed by Vazquez-Cendon [9] and Bermudez and Vazquez [10] for treating the homogeneous part of shallow water equations, which satisfies the exact conservation C-property. Integration of the source term over the control volume gives Following Bermudez et al. [11], this approximation is upwind, and the source term is replaced by a numerical source vector , such that Here, the point or is defined as the projection of the source term vector in the basis of eigenvectors of the Jacobian matrix. Thus the function source term is where is the identity matrix, is the Roe flux Jacobian, and represents an approximation of the source term at the cell interface .

To obtain higher-order spatial accuracy, the fluxes at each edge are calculated using a piecewise linear function of the state variable inside the control volume. For the cell-centred mesh, the MUSCL (Monotone Upstream Centred Scheme for Conserved Laws) approach is adopted [12], whereby the left and right values of the state variables are evaluated from in which is the vector distance between the barycenter coordinates of cells and and and represent the gradients on the triangles, respectively, and .

Time integration is undertaken using the second-order Runge-Kutta method [13] in which For stability, the global time step must satisfy the CFL stability condition given by where is the wave celerity. The Courant number CFL is set to 0.8 for all cases considered herein.

4. Simulations of the Nador Lagoon

Figure 1 shows the computational mesh fitted to the Nador Lagoon, based on interpolated measurements of the lagoon bed terrain. The mesh comprises a total of 9240 triangles and is refined according to boundary and bed gradient criteria, in order that the mesh density is more concentrated in the more complicated bed and near shore regions of the lagoon. Initially, the flow conditions are quiescent and the contaminant is distributed uniformly throughout the domain; that is,

At river flow inlet nodes, the contaminant concentration time history is specified. The water depth is set so that it is invariably positive, with mean value and fluctuating-free surface elevation such that where in which is the wave amplitude, is the angular frequency, and is the tide phase of the kth tidal constituent, noting that , , or in the present simulations.

For water exiting in the lagoon, the contaminant concentration is treated as having a transmissive boundary condition. For water entering the lagoon, the concentration at the open boundary nodes is set to zero, assuming that the incoming tide contains clean water. The water flow is forced with the main semidiurnal tide, of amplitude  cm which produces a maximum current speed of approximately 1.45 m/s through the lagoon inlet. Table 1 lists the physical and numerical parameters for this test case.

Table 1: Parameters and reference quantities considered in the present study.

The numerical model is used to assess the relative performance of two tidal inlet options in terms of certain key hydrodynamic characteristics of the Nador Lagoon, including water renewal and contamination dispersion. The existing inlet is denoted Inlet . The proposed new inlet is denoted . Wind effects are neglected. Uniform contamination at a concentration is imposed throughout the domain.

Figure 2 shows the time-dependent decay in average concentration over the entire lagoon, where the average concentration is defined by For both inlet options, about 70% of mass of the contaminant has been removed by the end of the sixth day. At all times, the average concentration obtained for Option remains consistently below that of Option . The results show that the tidal flushing is initially effective at decontaminating the lagoon but that a residual level of contamination persists after a week.

Figure 2: Time evolution of predicted average contaminant concentration .

Next, we define the normalized concentration as , where denotes the maximum value of the initial released concentration distribution . Table 2 presents the predicted results for water decontamination at three stations, owing to the tidal flushing through either Inlet or Inlet at three different times  h, 3 days, and 6 days after the initial contamination has occurred. Station 1 is located towards the north-west of the lagoon in a relatively stagnant subbasin, and it is not surprising that tidal flushing through either inlet has a much reduced effect. Stations 2 and 3 are located in the main basin, with Station 2 located nearer both Inlet and Inlet than Station 3. The normalized contaminant concentrations are quite similar between Stations 2 and 3 indicating that the tidal flushing at both inlets has a fairly uniform effect throughout the main body of the lagoon.

Table 2: Normalized contaminant concentrations obtained at three hypothetical sampling stations for the two inlet options O and N, at different times.

Figures 3 and 4 present predicted contamination concentration contours and depth-averaged velocity vectors for Inlet and Inlet , respectively. In both cases, the initial effect of decontamination in the lagoon due to the tidal forcing is evident in a zone near the inlet(s). Later, the decontamination occurs quite uniformly in the main basin, but zones of high concentration are evident in the northwest sublagoon region. Option seems better at reducing this contamination than Option . By the sixth day, the contamination has been considerably reduced throughout the basin, regardless of which inlet is considered. The tidal circulation patterns generated by Inlet comprise three primary counter-rotating gyres in the main basin, with secondary gyres in the northwest and southeast basins. For Inlet , the situation is similar, except that there are four counter rotating gyres in the main basin, and this tends to contribute rather more to mixing in the basin. In both cases, Station 1 is located near the middle of the northwest subbasin and close to a nearly stagnant zone within a gyre. Station 2 is near the shear layers delineating the boundary of one of the primary gyres in the main basin. Station 3 is in the shear layer between two counter rotating gyres, perhaps demarking a southwestern subbasin. The flow patterns are quite similar at each of the stations, as are the concentrations.

Figure 3: Spatial distributions of predicted (a) contaminant concentrations and (b) depth-averaged velocity vectors at three different simulation times (from top to bottom  h, 3 days, and 6 days) for tidal flushing through Inlet .
Figure 4: Spatial distributions of predicted (c) contaminant concentrations and (d) depth-averaged velocity vectors at three different simulation times (from top to bottom  h, 3 days, and 6 days) for tidal flushing through Inlet .

The numerical simulations indicate that Inlet is more effective at water renewal in the lagoon than Inlet . Obviously, the distribution of contaminants depends strongly on the tidal characteristics and the domain geometry. The present test cases demonstrate the effectiveness of the present finite volume model at representing steep horizontal gradients in the shallow flow and the plume like behavior of the initial decontamination zone within a complicated domain both in terms of bed terrain and boundary configuration.

The present finite volume solver required CPU time of 3 days and 11 hours to simulate tidal flows lasting for 6 days on a personal computer (whose configuration is described above).

5. Conclusions

A numerical model based on the conservative form of a coupled system of two-dimensional shallow water flow and advection-diffusion equations has been used to simulate contaminant transport in the Nador Lagoon, which has complicated bathymetry and an irregular coastal geometry. The simulations, though highly idealized, indicate the potential of the model in predicting contaminant transport at field scale and in evaluating the relative impacts on tidal flushing of two inlet options, and . It is found that situating inlet to the northwest of inlet leads to an additional gyre being created in the lagoon to the benefit of its decontamination by tidal flushing.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


  1. F. Benkhaldoun, S. Daoudi, I. Elmahi, and M. Seaïd, “Numerical modelling of sediment transport in the Nador Lagoon (Morocco),” Applied Numerical Mathematics, vol. 62, no. 12, pp. 1749–1766, 2012. View at Publisher · View at Google Scholar · View at Scopus
  2. F. Benkhaldoun, I. Elmahi, S. Sari, and M. Seaïd, “An unstructured finite volume method for coupled models of suspended sediment and bed-load transport in shallow water flows,” International Journal for Numerical Methods in Fluids, vol. 72, no. 9, pp. 967–993, 2013. View at Publisher · View at Google Scholar · View at Scopus
  3. A. Amahmouj, E. M. Chaabelasri, and N. Salhi, “Computations of pollutant dispersion in coastal waters of Tangier's Bay,” International Review on Modelling and Simulations, vol. 5, no. 4, pp. 1588–1595, 2012. View at Google Scholar · View at Scopus
  4. F. Boushaba, E. M. Chaabelasri, N. Salhi, I. Elmahi, F. Benkhaldoun, and A. G. L. Borthwick, “A comparative study of finite volume and finite element on some transcritical free surface flow problems,” International Journal of Computational Methods, vol. 5, no. 3, pp. 413–431, 2008. View at Publisher · View at Google Scholar · View at Scopus
  5. R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, Cambridge, UK, 2002.
  6. P. L. Roe, “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational Physics, vol. 43, no. 2, pp. 357–372, 1981. View at Google Scholar · View at Scopus
  7. J. G. Zhou, “Velocity-depth coupling in shallow-water flows,” Journal of Hydraulic Engineering, vol. 121, no. 10, pp. 717–724, 1995. View at Google Scholar · View at Scopus
  8. F. Alcrudo, P. Garcia-Navarro, and J.-M. Saviron, “Flux difference splitting for 1D open channel flow equations,” International Journal for Numerical Methods in Fluids, vol. 14, no. 9, pp. 1009–1018, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. M. E. Vazquez-Cendon, “Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry,” Journal of Computational Physics, vol. 148, no. 2, pp. 497–526, 1999. View at Publisher · View at Google Scholar · View at Scopus
  10. A. Bermudez and M. E. Vazquez, “Upwind methods for hyperbolic conservation laws with source terms,” Computers and Fluids, vol. 23, no. 8, pp. 1049–1071, 1994. View at Google Scholar · View at Scopus
  11. A. Bermudez, A. Dervieux, J. A. Desideri, and M. E. Vazquez, “Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes,” Computer Methods in Applied Mechanics and Engineering, vol. 155, no. 1-2, pp. 49–72, 1998. View at Google Scholar · View at Scopus
  12. I. Elmahi, Schémas volumes finis pour la simulation numérique de problème à fronts raides en maillages non structurés adaptatifs [Ph.D. thesis], Université de Rouen, Mont-Saint-Aignan, France, 1999.
  13. Ch. Hirsch, Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics, vol. 1, John Wiley & Sons, New York, NY, USA, 2nd edition, 2007.