1. Introduction

Lagoons are almost closed bodies connected loosely to the sea by narrow channels subjected to siltation and sediment obstruction. The hydrodynamic circulation pattern of these bodies is difficult to compute due to their usually complex geometry combined with the influence of the tides, the wind and the fresh water inputs coming from the rivers and streams of the lagoon watershed. From a biological point of view, lagoons can be characterized as complex ecosystems in a fragile equilibrium and their modeling is a difficult task as one has to take into account numerous interactions between physical, chemical and biological parameters.

Lagoons have important environmental and societal aspects. They are extremely rich from a biological point of view as the reproduction area of a great number of species as well as the passage of many migratory bird species. Lagoons are also the siege of many human activities: urbanization, tourism, agriculture, fishing and shell-fishing farming all influencing the lagoon state in a complex way. Many lagoons are affected over the years by anthropogenic influences, including environmental degradation and local or global sources of pollution.

The modelling of lagoon poses a lot of interesting and difficult questions to numerical analysis: The geometrical complexity of these areas requires the use of sophisticated two or three dimensional finite volume or finite element methods that are computationally very expensive and ask for adaptive meshes and parallel computing. Beside this, lagoon’s modelling is intrinsically a multi-physics problem. It therefore involves multiple spatial and temporal scales that have to be connected through parameterization as well as by well designed numerical tools.

Among these lagoons, the Nador lagoon located on the North-East of Morocco is an ecosystem that is of great biological interest, ecological and economic. It covers an area that exceeds 120 Km² and a maximum depth of 8m (see Figure 1), and is fed by water of the Mediterranean through a pass known as ‘Bokhana’, the freshwater waterways, the rejections of the untreated human activities (agriculture and urban water industry : metallurgy, textile …), and by water of the station of purification.

Figure 1. Schematic description of the Nador lagoon

The Nador lagoon is among the most productive areas of the Mediterranean. It has a particular interest in the region, both ecologically, economically that and exploitation. It also presents significant opportunities for economic development and tourism in the region, thanks to its geomorphological structure characterized by its easy access, its proximity to seaside resorts and its proximity to urban areas (Idrissi et al. [1]). The lagoon has attracted for decades, the attention of many scientific organizations and investments (INRH, universities, etc). Its rich variety of flora and fauna makes ecological environments play a key role in the environmental balance. However, ecological qualities of the lagoon have been relatively affected, during the last years, by the human activities, particularly the pollution coming from worn water of the city and the waste evacuated along its banks as well as the hydrocarbons coming from the industrial boats.

In addition, the Nador lagoon has been the subject of many recent investigations that have concerned water quality, currents, bathymetry, flora, fauna, fishing and aquaculture. In 1981, Tesson and Gensous [2] established the first distribution map of the clay fraction at wide lagoon. Thereafter, using statistical tools, they studied the mode trapping trace elements in different mineral fractions. In 1987, Irzi [3] studied the sedimentological and micro-faunistic aspect of the lagoon. She highlighted the problem opening-closing of the pass and the distribution of foraminifera and ostracods depending on intra-lagoon currents and anthropogenic factors. In 1995, Inani [4] established a diagnosis of the state of pollution in sediments and waters of the lagoon. He highlighted four issues: anoxia, sursalure, accumulation of the organic matter and, finally, high metal concentrations in micro-pollutant. In 1996, Lefebvre and his co-workers [5] published a study which focused on the bio-geological evolution of the Nador lagoon during the period of 1982-1993. They concluded that the main change in the ecosystem is eutrophication and the significant increase of organic matter by urbanization and aquaculture activities. These conditions contributed to a one-time silting of the ecosystem that has resulted in a reduction in the extent of the deposit of the clam. In 2002, Irzi [6] deepened his studies on the ecology of the lagoon. She showed that the distribution of foraminifera is low at the continental margin (contribution of urban waste) and near the fish farms (very rich areas of organic matter). The same year are presented studies on heavy metal contamination in intra-lagoon sediments and especially in the vicinity of the city of Nador [7]. In 2005, Bloundi [8] conducted a geochemical study of the Nador lagoon including the identification of anthropogenic impacts on water and lagoon’s sediments. This study has allowed estimation the heavy metal contamination and eutrophication of the lagoon ecosystem of Nador. In 2011, Kada [9] realized a distribution map of the salinity of surface water in the lagoon of Nador, followed by biological and ecological study and characterization of populations in the lagoon. He also showed that the salinity of the lagoon remains close to the salinity of the Mediterranean Sea, and that this salinity is directly related to metrological conditions, continental inputs, as well as communication with the sea via the pass.

As one can see, almost all the studies realized on the Nador lagoon have focused especially on the environmental aspect of the lagoon and are biological, geochemical and economic. However, the studies from numerical point of view concerning the hydrodynamic, water circulations and transport-dispersion of contaminants in the surroundings of the lagoon are missing.

In this paper, we propose a finite volume method for the simulation of the hydrodynamic and water circulations in the Nador lagoon. The physical model is based on the 2D shallow water equations that include slope variations and friction terms. The computational code uses a Non Homogeneous Riemann Solver (SRNH) for the discretization of the convection and bed slope source terms [10], and semi-implicit splitting scheme for the friction terms [11]. It is of second order in space and time and can be used on complicated geometries by incorporating unstructured meshes. In addition, the algorithm preserves the balance between the pressure gradient and source term due to the slope variation. The preliminary results that will be presented here considering the new “pass” with the Mediterranean Sea, will demonstrate the performance and robustness of the solver for the simulation of flows in real situations.

The structure of this paper is as follows: section 2 presents the set of Saint-Venant governing equations, section 3 describes the basic principles of the numerical methods used, the extension to second order in space and time and the treatment of the friction terms are also shown. In section 4, we present an application of the solver to the simulation of the hydrodynamic circulations in the Nador lagoon. Some concluding remarks are given in section 5.

2. The Equations to be Solved

The mathematical model can be derived by vertical integration of the three-dimensional incompressible Navier-Stokes equations along with the assumptions of a hydrostatic pressure and a vertically uniform horizontal velocity field, which results in the well-established shallow water equations written in conservative form as

(1)

W is the vector of conserved variables, S₁ and S₂ are the source terms due to slope variations and friction forces, and F and G are the advection tensor fluxes.

, ,

,

where x and y are the spatial coordinates, t is time, h is the height of the flow above the bottom Z, g is gravity acceleration, u and v are the velocity components in the x and y directions, and the water density. and are the bed shear stress in the x and y direction, respectively, they are defined by the depth-averaged velocities as

where C_b is the bed friction coefficient, which may be either constant or estimated as , where is the Chezy constant, with being the Manning coefficient at the bed.

3. Numerical Method

3.1. Mesh Generation and Bed Topography

Let Ω be the domain of computation covered by a conformal grid consisting of a set of unstructured polygons K that are as general as possible. In this work the unstructured meshes are composed of triangles. Figure 2 shows the mesh of the Nador lagoon generated by using the Delaunay triangulation technique [12]. Since the boundary conditions in the Mediterranean part are not well known, the computational domain has been limited in this study to the pass between the lagoon and Mediterranean sea. In this way the tidal boundary conditions will be imposed at this pass.

Figure 2. Unstructured mesh of the Nador lagoon

The geometry and the bed surface topography of the lagoon are very irregular and several regions of various depths coexist. In our simulations the bathymetry was reconstructed from topographical data. This bathymetry is illustrated in Figure 3 below.

Figure 3. Bathymetry of the Nador lagoon

3.2. Finite Volume Formulation

Research on numerical solution of equations (1) has received considerable attention during the last decades and a several finite volume methods have been developed, compare [13-15] among others. The main advantages of these methods lie on their implementation on unstructured triangular meshes and preserving conservation properties of the equations. The computational code developed here uses the «cell-centred» finite volume formulation for which all the state variables are updated at the centroid of each cell (the centre of gravity in our case). Figure 4 illustrates the type of control volumes used.

Figure 4. A cell-centred control volume

Hence, integration of the system (1) over a control volume T_i and using the Gauss divergence formula, a finite volume discretization yields to

(2)

where is the convective flux tensor defined by

N(i) is the set of neighbouring triangles of the cell is the edge separating control volumes T_i and T_j, and n=(n_x,n_y) denotes the unit vector normal to

3.3. Discretization of the Convective and Bed Slope Source Terms

The system of Saint-Venant equations (1) includes the variation of topography and therefore special attention should be paid to his numerical discretization. Indeed, a well known problem is that shallow-water equations on non-flat topography have steady-state solutions in which the flux gradients are non zero but are exactly balanced by the source terms. Standard numerical methods for the discretization of conservation laws may fail in correctly reproducing this balance (C-property, see e.g., [16]) and, thus, specific methods have been developed to deal with this problem (well-balanced schemes, see e.g., [16-18]).

The method adopted herein for the space discretization has been proposed by Benkhaldoun, Elmahi, and Seaid in a series of papers (see for instance [10, 19]). It is based on a Riemann Solver, named SRNH, and consists of a predictor stage for which the state W_ij at each edge is computed by solving a Riemann problem, and a corrector stage where the state Wⁿ⁺¹ at time tⁿ⁺¹ is reconstructed iteratively using the physical flux on W_ij.

The scheme considers only the advection part and the source term containing slope variations

(3)

Let's start with the predictor stage. A projection of the system (3) according to the normal η and tangential yields:

(4)

where and are, respectively, the normal and tangential velocity.

For this projected system, the predictor stage is then formulated, for the evaluation of the averaged state on each edge , using an upwind scheme, in the following manner:

(5)

where

, ,

In (5), is the intermediate Roe-averaged state given by

denotes the sign matrix of the Jacobean, it is defined by:

and

and are respectively the eigenvector and eigenvalue matrices of .

By incorporating these matrices in the predictor stage (5), the projected state on each edge can be easily obtained. The conservative state is then evaluated using the transformations and.

Therefore, the incremental corrector stage is written using the non projected conservative states in the following manner:

(6)

being the surface of the control volume and the length of the edge separating triangle and .

In the corrector stage (6), the source term approximation is reconstructed in such a way to obtain a well balanced scheme preserving steady state at rest (see [10] for details on this reconstruction).

3.4. Second Order Extension

The SRNH scheme described above is first order accurate in space. It is monotone but has a poor accuracy due to the large amount of numerical dissipation. The extension to 2nd-order accuracy in space can be achieved by using a classical MUSCL technique [20]. In the definition of the flux in (6), we replace the piecewise constant values and by more accurate reconstructions deduced from piecewise linear approximations, namely the values and , reconstructed on both sides of the interface as follow/

(7)

where denotes the gradient operator, and are respectively the barycentres of cells and . is a parameter between 0 and 1. In practice one uses The resulting scheme is second order in space but is not necessarily monotonous and non-physical oscillations are produced. To damp the numerical oscillations in the current computations, the Van Albada flux limiter is applied:

(8)

with the Van albada limiter given by

(9)

where . With this limitation one obtains a monotonic second reconstruction for the hyperbolic part.

3.5. Treatment of the Friction Terms

The friction terms in (2) are discretized using an operator splitting procedure. For example, to evaluate the x-momentum, the second equation of the system (2) is split into two equations

(10)

where is the Manning’s coefficient and Res describes the convection contributions in the x-momentum equation corresponding to the surface integral in eq. (2), and is approximated as the sum taken over all boundary segments.

First, a semi-implicit method is used to integrate the upper equation in (10), giving

(11)

In the second step, the value is taken to be the initial condition when solving the second equation in (10).

3.6. Explicit Time Advancing

The spatial discretization of the convection and source terms leads generally to the following semi-discrete form:

(12)

where is an approximation of the following discrete contributions:

To obtain an explicit scheme first-order accurate in time, the time-discretization can be carried out using an explicit Euler method, so that the resulting first-order numerical scheme is

(13)

Similarly, in order to obtain an explicit scheme second-order accurate in time, the time discretization can be carried out using a second-order TVD Runge–Kutta method [21]. The expression of the second-order numerical scheme is

(14)

4. Numerical Test

The test case that we consider here consists in the study of the hydrodynamic and water circulations in the Nador lagoon. This test problem is interesting in the sense that it can give an answer to how the water moves into the lagoon, regions where the water recirculates and those where it stagnates. These last regions can be in fact very affected by contamination of wastewater coming from the city. The computations can also give an idea of the time residence of the water coming from the Mediterranean Sea before leaving the entrance of the lagoon.

Initially, the flow is assumed to be at rest with constant free surface (h + Z = Cte).

Two types of boundary conditions must be specified respectively at the pass between the Mediterranean Sea and the lagoon noted and land-type boundary located on the Nador coastlines Hence the following boundary conditions

on and

are imposed for the water flow and

(15)

are used for the water height.

In (15), H is the depth from a fixed reference level to the bottom, h₀ is a given averaged water elevation taken here equal to 3m (initially we have h=H+ h₀), is the tidal amplitude at the entrance of the lagoon, is the angular frequency of the tide, and is the phase of the tide. In the Nador lagoon, the main astronomical tidal constituents are the semidiurnal M₂, S₂ and N₂ tides, and the diurnal K₁ tide. The values of the parameters for each tidal wave are taken from [22] and are illustrated in table 1. Thus, the boundary condition on is replaced by

(16)

Note that the boundary conditions for the water height on are time-dependent and should be updated at each time step according to the condition (16).

Table 1. Parameters for the considered tidal waves at the “pass” of the Nador Lagoon

The computations have been performed with a physical time step chosen in such a way that the following stability condition is satisfied:

where is the eigenvalue of the system, evaluated at the interface between two cells T_i and T_j, and CFL is the courant number taken here equal to 0.5 in order to ensure stability of the numerical scheme.

In the figure 5, is presented the water height inside the lagoon. The computations have been performed for 6 days in physical times.

Figure 5. Water depth contours in the Nador lagoon

In Figures 6 and 7 respectively, are illustrated the velocity field together with the isolignes of u-velocity in the case of high and low tides. The clear result is the presence of recirculation zones inside the lagoon. It is also clear that in case of high tides the water flow is coming from the Mediterranean part and feeds the lagoon, this can be clearly seen on the Figure 6. On the other hand, when the tides are low, the water flows from the lagoon toward the Mediterranean part (see Figure 7). One can also observe that the flow is almost inert at Beni Enzar and Kariat Arkmane regions. These regions could then be very affected in case of contamination by wastewater because the pollutant can stagnate in these regions.

Figure 6. Velocity field (top) and isolignes of u-velocity (bottom) in the case of high tides

Figure 7. Velocity field (top) and isolignes of u-velocity (bottom) in the case of low tides

In order to validate the results obtained by our solver, and since the observation data are missing, computations have been compared with those obtained by using the finite volume Vazquez solver which is based on a well balanced modified Roe scheme (see [16]). Figures 8, 9 and 10 show the historic of water depth and u-velocity with time at three gauged points, G₁ situated at the entrance near the “pass”, G₂ situated inside the lagoon and G₃ situated at Kariat Arkmane region where the flow is almost inert (see Figure 5). One remarks that the velocity is larger at the gauged point G₁. This is due to the tide coming from the pass. The values of this velocity evolve from negative to positive depending on the case of high or low tides. Moreover, the negative and positive values of u-velocity at the gauged point G₂ can be explained by the hydrodynamic circulations in this zone. Finally, at the gauged point G₃ which is far from the tides and recirculation zones, we remark that the velocity is almost zero. One can notice also an excellent agreement between our results and those obtained by the well known Vazquez finite volume solver, which confirms the capability of the finite volume solver developed here to deal with complicated geometries and real situations in shallow water flows.

Figure 8. Water depth and u-velocity versus time at gauged point G1

Figure 9. Water depth and u-velocity versus time at gauged point G2

Figure 10. Water depth and u-velocity versus time at gauged point G3

5. Conclusions

In this paper, a numerical code based on the finite volume method has been applied for the simulation of the hydrodynamic and water circulations in the Nador lagoon. The method consists of two stages, which can be interpreted as a predictor–corrector procedure. In the first stage, the scheme uses the projected system of the coupled equations and introduces the sign matrix of the flux Jacobian, which results in an upwind discretization of the characteristic variables. In the second stage, the solution is updated using the conservative form of the equations and a special treatment of the bed bottom to obtain a well-balanced discretization of the flux gradients and the source terms. The solver is fully second order in space and time and can be used on complicated geometries with unstructured meshes. The test case presented here shows the good performance of the method and confirms its capability to provide accurate and efficient simulations for shallow water flows including complex topography and friction forces on unstructured grids.

The next step of this work will consist in the study of the transport-dispersion of contaminant produced by the wastewater in the Nador lagoon. The results obtained in this paper for the hydrodynamic will be considered as an initial condition for the simulation of pollutant transport problem. A relevant study would also measure the residence time of water coming from the Mediterranean Sea into the lagoon in order to find a way to renew often, and thus make it less polluted. An adequate numerical study would determine the necessity and indeed the eventual location of another passes between the lagoon and the Mediterranean permitting to reduce the residence time of a given tracer.