Abstract
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.
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.
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.
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.
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.
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.