Two-Dimensional Numerical Modeling of Morphodynamic Evolution in Bouregreg Estuary (Morocco)

The Bouregreg valley, located on the Atlantic coast of Morocco, has been subject to several anthropogenic and natural impacts that disturb the proper functioning of the estuary’s ecosystem. Important dredging operations influ - enced by both tides and freshwater inflows have created a significant variation of estuary’s morphology character - ized by a sand bar developed at the mouth. The main objective of this study is to investigate the Morphodynamic processes of the estuary and in particular the impact of dredging operations using a two-dimensional horizontal numerical (2DH) modeling approach based on the shallow water equations. Several numerical tests have been carried out to calibrate the model and to investigate the hydrodynamic and morphological aspects of the estuary. For this technical solution, the modelling simulations show that in addition to the quantities of sand that will be dredged initially, the maintenance of the channel will facilitate a better water circulation between the Bouregreg valley and the Atlantic Ocean and improve the stability of the civil engineering structures and the neighbouring monuments. In addition, this water circulation will guide the river flow in the ebb phase of the tide. The model results show a reduction of about 30% in the impact of erosion, observed before the digging of the cana l.


INTRODUCTION
The Bouregreg valley area (Atlantic coast of Morocco) located between Sidi Mohamed Ben Abdellah dam and the mouth of the estuary, is a complex system characterized by several physical and environmental aspects and a diversity of flora and fauna. In the past, the Bouregreg estuary has known dynamic navigation activities related to trade and fishing. The old port of Rabat-Sale was at the beginning of the 20 th century, the most important commercial port in Morocco. The decline in traffic of this port due to the port of Casablanca extension and the importance of dredging costs for the maintenance of its channel, led to the cessation of its activity in 1940 and since then, only the activities of pleasure and artisanal fishing have been maintained. Major morphological modifications have been caused after the construction of the Sidi Mohamed Ben Abdallah dam in 1974 located approximately 23 km from the mouth of the estuary, especially hydraulic and sediment processes. The large quantity of sediment particles transported by the Bouregreg river, which constituted the main supply resource in terms of sediment budget, is stopped by the dam upstream of the river as well as an impact on the frequency and intensity of flowing and Two-Dimensional Numerical Modeling of Morphodynamic Evolution in Bouregreg Estuary (Morocco) the aquatic habitat of the Bouregreg river. This area is subject to multiple actions of environment degradation: wastewater discharges into the river, solid waste deposits, air pollution. Many aspects due to the proximity with Rabat and Sale cities, the valley and particularly its downstream part, present high traffic exchange, and huge urbanization and industry development. On the hydraulic, morphological and sedimentological levels, the river reach is subject to periodic effect of Atlantic Ocean's tides and occasional flows due to the dam water releases. At the estuary mouth, the penetration of ocean waves contributes to the formation of a sandy bar and the estuary is constantly subject to the influence of the tides that reach up to the dam. This situation threatens the stability of important infrastructures, and disturb fish activities. Therefore, in order to identify all these phenomena impacts and to anticipate and control their repercussions, Moroccan's authorities develop a management plan for the valley of Bouregreg. For this purpose, different specific technical environmental and economic studies have been carried out to better understand these processes and improve the situation.
This research paper is a contribution in that way to evaluate the hydrodynamical and sedimentary aspects of the Bouregreg estuary to study the impact of sediment extraction on the morphodynamic of the estuary. This approach is based on a global two-dimensional numerical approach coupling hydrodynamics, sediment transport and morphologic movements. The main objectives are to assess the response of the estuary's bathymetry to the excessive dredging of sand; and to help managers and decision-making to enhance technical and environmental solutions to improve the situation. Many numerical simulations are used to evaluate the impact of sand dredging on the spatial and temporal evolution of the estuary's bathymetry and to analyze a compensatory measure representing the dredged channel solution and its efficiency to improve the bathymetric situation.

THE BOUREGREG ESTUARY Overview
The Bouregreg river, is a part of Bouregreg watershed located on the Atlantic Ocean at 34°North and 6°50W and covering an area of 9700 km 2 . The Bouregreg estuary on the edge of the cities of Rabat and Salé, with a length of 24 km long and an average width of 150 m, is artificially limited upstream by the Sidi Mohamed Ben Abdallah dam, to supply drinking water of many cities of the Moroccan Atlantic coast ( Figure 1).

Geological aspect
The Bouregreg estuary is located in the coastal Meseta characterized by plains and plateaus where the altitude decreases as one approaches the Atlantic Ocean. Its drains, in its upstream part, a substratum constituted of Paleozoic lands of Schistosandstone and Dolonutic nature in which rocks of magmatic or volcano-sedimentary origin can sometimes be intercalated (Beaudet, 1969).

Climatic aspect
The climate of the Rabat region is Mediterranean, characterized by a dry season from May to October (average over 25 years: 0.1 mm in July), hot summer and Autumn with high temperatures varying between 18°C to 27°C, and a wet period relatively cold varying between 8°C and 16°C. We observe an irregular distribution of rainfall varying between 370 mm and 800 mm per year characterized by variability and brutality of rainfall. The marine winds (SSW to W), are predominant and regular (3m/s) (CID et al., 2003;Sadok, Charafi, Kamal, 2003).

Hydrological aspect
The Bouregreg river with a length of about 300 km has an annual flow of about 7.10 m 3 /s estimated before the construction of Sidi Mohammed ben Abdellah dam. In the actual situation, the flow became negligible. The hydrodynamics of the estuary is essentially dependent on the semi-diurnal tidal regime. The saline intrusion is characterized by a decreasing gradient from the dam depending on the intumescence. The Bouregreg river flow is very low compared to the tidal currents, so the dispersion becomes responsible for most of the exchanges within the estuary. Therefore, the estuary is in major time with homogeneous salinity moving gradually to the upstream due to the freshwater water inputs decrease in time especially from Sidi Mohammed ben Abdellah dam (CID et al., 2003;Sadok, Charafi, Kamal, 2003) Past studies of Bouregreg estuary Since many decades, the Bouregreg estuary has been the subject of several academic researches, consultancy studies and institutional reports. These works are generally based on empirical / semi-empirical approaches or simplified analytical and numerical approaches investigating different processes of the Bouregreg estuary. We can cite for the water quality aspects: intrusion salt, contamination of sediments by heavy metals, impact of different kind of pollutions on the water quality (biology, chemistry, biogeochemistry). For hydrodynamic and sediment studies, there is few works using simple approaches investigating for example the tide flow, water circulation. All these information's and outcomes are important to better understand the functioning of the estuary's ecosystem especially after the construction of Sidi Mohammed ben Abdallah dam. However, there is a luck of studies based on numerical approaches applied to the Morocco's estuaries in general and particularly for this area. The authors of this paper developed in the past two decades a global process-based numerical model HYSED in both two-dimensional horizontal and quasi-threedimensional ways which can be adapted for different situations and water systems (lakes, estuaries, and river) depending on the main objectives of the study considered, the type of sediments and the quantity and quality of database available. For sediment and morphological processes, it can investigate both suspended sediments resolving the advection-diffusion of sediments equation and its variation with depth of the flow with a real time coupling the two phases of the mixture, or total sediment transport and the morphology of the bed tacking into account the sediment's cohesion

MODELING OF THE BOUREGREG VALLEY'S HYDRODYNAMICS AND MORPHOLOGICAL PROCESSES
Numerical models are important tools that help to study estuary's problems, allowing managers and decision-makers to simulate different conditions, and to predict the evolution of estuaries in the future and to choose and optimize management solutions. These models applied to estuary's systems have known important development due to the applied mathematics and computer sciences huge progression in last decades. They become important tools to better understand the impact of different natural and artificial aspects: physical, ecological and socio-economic which are involved in these systems. Coupled to in-situ measurements and data bases development, this hybrid approach become unavoidable to better design and implement management policies and to optimize decision-makers choices for sustainable development of these complex ecosystems  The modular structure of the numerical model is based on a two-dimensional finite difference approach, coupling the Hydrodynamic module (HY_2D) and a sedimental module (SED_2D). This real time coupling between the two modules in each node of the grid and at any time of simulation, allowed the numerical model to simulate the evolution of the bottom and its configuration. The model provides at each node of the mesh domain and after each time step: the velocity field of the fluid flow depth-integrated in x and y directions (U and V), the water depth (H), total sediment loads, and the variation of the bottom bathymetry (Z b ). The interaction of the hydrodynamic flow with the bed level variation is taken into account in the next step of modeling using a dynamic ap-proach ( Figure 2).

Description of the hydrodynamic module
The hydrodynamic module (HY_2D) is based on the two-dimensional horizontal numerical resolution of the 2D shallow water equations depthintegrated (generalized Saint-Venant) (  (1) where: H -water's depth, U and V -respectively the components of flow velocity depth-integrated in x and y directions, t -time, g -acceleration of gravity, f -Coriolis parameter, Z b -bottom level, K -an aerodynamic coefficient, W -wind speed, ϕ -angle between wind speed vector and x axis, n -manning coefficient, ε -flow turbulent viscosity coefficient.
The Finite difference equations are formulated according to an improved version of the Mac Cormack fractional scheme (Garcia-Navarro and Saviron, 1992; Paola and Voller, 2005). This scheme is based on an explicit fractional step method where the finite difference operator is decomposed into a sequence of simpler operators. In order to numerically integrate the equations, the domain under consideration is superimposed with a grid where all the variables are defined at the centers of the cells in order to represent the mean properties of each cell.

Description of the sediment transport module
The sediment transport module (SED_2D) is composed of two sub-modules: CTS_2D and F_2D) and simulates the variation of the bottom level and the spatial and temporal evolution of the valley's morphological (erosion, transport, deposition).
CTS_2D provide the components of the total sediment transport load in the directions x and y. It is based on different total sediment transport equations which are recommended by the specialized literature, especially Engelund-Hansen, Van Rijn and Ackers -White (Rijn, 1987;Sadok, 1988; Van der Wegen and Roelvink, 2008;Van Rijn, 1993). For the current case study, after analyzing the hydraulic and granular conditions and tacking into account the luck of some specific data of Bouregreg estuary, we choose the Engelund-Hansen formula (Engelund and Hansen, 1967;Sadok, 1988). This formula which include implicitly bed load and suspension load is characterized by its simplicity and accuracy as tested and recommended by different authors. However, the procedure of 2D version of this sediment formulae was used to take into account the slopes of the bathymetry in each cell of the grid (Van der Wegen and Roelvink, 2008; Van Rijn, 1993).
The Engelund-Hansen formula (Eq. 4a) based on an energy conservation concept and relates the total sediment transport directly to the water velocity and the median gran size of sediment particles is: where: q t -magnitude of the total load transport (m 3 /ms); q s -magnitude of the suspended load transport (m 3 /ms); q b -magnitude of the bed load transport (m 3 /ms); U m -magnitude of the integrated velocity over the vertical direction (m/s); C -Chezy friction parameter (m 1/2 /s) (Eq. 4b); s -sediment density, default 2.65); d 50 -median grain size of bottom sediment particles (m); S -bed slope; H -water depth (m).
F_2D is a morphological sub-module which provide the evolution of the bed of the flow. It is based on a solid mass continuity equation (generalized Exner's equation) (Eq. 5) representing the masse conservation of the sediment mixture inside an elementary control volume and a balance between the divergence of total sediment transport field and the evolution of the bed level corrected for bed porosity will allow to write Sadok and Marche, 1988; Van der Wegen and Roelvink, 2008;Wang et al., 1995).
where: q t,x and q t,y -the components of the total sediment transport per width along x and y, Z b = the mean bathymetric level, p = the porosity, default 0.4.
The numerical model simulates the temporal evolution of the mean bathymetric level of the bottom Z b , at each node of the mesh, representing the study area. To resolve numerically this mass conservation equation of the sediment bed mixture we choose the Lax numerical scheme in order to reduce the numerical diffusion and obtain a stable solution. This approach was used in many studied with good satisfaction (Van der Wegen and Roelvink, 2008; Van Rijn, 1993).
The equation (Eq. 6) discretized by means of the diagram of Lax is written: where: Z i,j -bottom elevation at each node (m); Δt -time step; Δx and Δy -space step; λnumerical smoothing parameter.
The parameter λ determines to what extent the bathymetric levels surrounding the level Z i,j at time (t) are taken into account for the calculation of a new bottom elevation Z j,i at time (t+Δt). This allows for numerical smoothing at sharp transi-tions in the bathymetric surface (Van der Wegen and Roelvink, 2008).
The boundary conditions represent the sediment transport discharges at open and closed boundaries. At open boundaries, the sediment transport discharge is given as a time series of data (sediment rating curve) or as function of time q t = f(t). At closed boundaries, the condition of zero normal (boundary) flow is considered.

Description of the study area
The study area represents the downstream estuary part, from the mouth (ocean's inlet) to the first bridge encountered on Bouregreg river, with a length of about 2 km (Figure 1). Because of difficulties to define boundary conditions in the bridge section of the river the mathematical model mesh was extended over a length of about 23 km, to reach the dam's structure where both the tide and water releases, are measured.
The estuary's area studied is represented by a schematic configuration (discretized grid) of 1500 meshes with 2500 m 2 each one (50X50). The nature of the tide flow generated imposes a wetted area variable in time, depending on the phases of the tide. To solve this problem, an automation of the calculation which allows to alternate two configurations related to the tide's flow and tide's ebbing phases has been adopted. Thus, the number of calculation points (nodes) on the mobile grid follows the evolution of the coastline mooring and reaches respectively 1500 and 1160 points during the ebb and flow phases. The length of the domain starts from downstream to upstream as shown in the plan view of the area (Figure 3).
The numerical model developed takes into account the water releases from Sidi Mohammed ben Abdellah dam (upstream end) and the variation of the water level at the mouth of the Bouregreg estuary (downstream end). The time step used is 10 seconds and the model simulates at each time step, the variation of the water level, the average velocities along the main directions (X and Y) and the resulting flows. It also calculates the maximum stresses exerted on the bottom. Also, the model defines the maximum currents, reached respectively during the flow and the ebb phases in each node of calculation. Finally, the model allows a readjustment of the ratings at each mesh and gives consequently the spatial and temporal evolution of the bathymetric state of the estuary.

The data base
For modeling the Bouregreg hydrodynamics and morphologic processes we used multiple available data especially: bathymetry of the area, granulometry of sediment particles, hydrodynamics of the estuary (tides and currents), water releases from Sidi Mohammed ben Abdellah dam.

Bathymetry
The initial topographic and bathymetric configuration were derived from a campaign of measurements realized over the area separating the mouth of the estuary and the bridge (Sadok, Charafi, Kamal, 2003). The results are shown in (Figure 4).

Sediment granulometry
Using many sampling technics for granulometric and sedimentological analysis to investigate the granulometry of sediment particles of the concerned part of the estuary allowed us to choose an average diameter d 50 of 0.25 mm as representing sediment particles (Sadok, Charafi, Kamal, 2003).

Hydrodynamic of tides
The tide in the Atlantic ocean's inlet (mouth) is semi-diurnal with an average period of approximately 12h25mn; an average amplitude of 2.80 m at spring at the ocean's inlet and a range varying from 3.80 m to 0.5 m at spring and neap respectively. The velocity of tidal propagation along the Bouregreg river decreases from downstream to upstream with a sinusoidal shape. The amplitude of the tide is modulated with a periodicity of 14 days corresponding to half -moon season. The hydrodynamic characteristic of the tide is listed in the Table 1. The tidal range is maximal at spring tide and minimal at the neap tide.
As a first approximation, the tide is approximated in this study ( Figure 5) at the level of the two bullheads (Eq. 7), by: where: H(t): water level /ZH; t -time; T 1 -tides' period (12h 25 min); T 2 -modulation period (14 days).

Dam's water releases
The drought experience in the last decade in Morocco and the high and growing demand for drinking water, and the luck of other feasible and economical solution, push water managers to raise the height of sidi Mohammed ben Abdellah dam because it represents a strategic surface water resources infrastructure for supplying many cities on the Atlantic coast. Therefore, the frequency and magnitude of releases of water from the dam will decrease in future especially due to droughts and climate change impacts. However, it is important to examine the influence of the most representative releases on the dynamics of the estuary. The monthly releases beginning at the 8th day, commonly made are shown in (Figure 6).

Sensibility analysis of relevant parameters
The hydrodynamic module (HY-2D) of the numerical model HYSED involves two fundamental parameters that are difficult to control and to evaluate analytically. In on hand the Manning's coefficient (n) tacking into account both roughness due to grains and to sand waves (Sadok, Charafi, Kamal, 2003;Sadok, 1988;Sadok and Marche, 1988;Van Rijn, 1993); and in the other hand the coefficient of turbulent diffusion (ε). A sensitivity study on ε showed that the maximum error generated for a variation of this coefficient from 0.1 to 1 is 6 to 12%, between measured and simulated values. The Table 2 gives the maximum water velocities obtained respectively during the flow and ebb phases as well as the variations generated by the various values of parameters for a constant configuration of the bottom.
However, the Manning's coefficient variation has clearly influenced the nature of the flow. Indeed, a variation of ''n'' from 0.02 to 0.04 generates differences that can reach 40 to 60%. In the Table 3 are listed the maximum flow velocities reached respectively during the flow and ebb phases, as well as the differences generated for different values of n.

Hydrodynamics' modeling
The global numerical model is calibrated to well reproduce hydrodynamic and morphology of the estuary. As shown by the sensitivity analysis, the Manning coefficient n is a very sensitive parameter, justifying its choice for the calibration of the hydrodynamic module. The Table 4 highlights the maximum of water velocities simulated by the model for different values of the Manning parameter, respectively during the flow and ebb of the tide. This Table 4 shows also the maximum   (Sadok, Charafi, and Kamal, 2003).
The comparison between the maximum velocities simulated by the model and those measured, allowed to adjust the model with a Manning's coefficient of 0.029. Figures 7 and 8 show the water circulation velocities simulated during the filling and emptying of the estuary. Figure 9 shows a significant phase shift and an attenuation of the tide's amplitudes from the mouth of the estuary which decrease when we move away from the Atlantic's Ocean. Therefore, currents induced in this area remain important, reaching 0.92 m/s and 1.16 m/s respectively during the ebb and flow phases of tide. These currents have the capacity to influence the hydrodynamic and sedimentary processes of the estuary and to impact consequently the equilibrium of the estuary's bottom bathymetry.

Numerical modeling of the impact of estuary's dredging
The purpose of the simulations carried out by the numerical discussion model is to evaluate the negative impacts of dredging operations on the hydrodynamic and sedimentary dynamics of the estuary. Thus, we simulate the spatial and temporal evolution of the estuary's bottom under the simultaneous action of tide's propagation and water releases from sidi Mohammed ben Abdellah dam. Figure 4 represent the initial estuary's bathymetry  of the estuary. Figures 10 and 11 shows the bathymetry's configuration simulated by the mathematical model and its variation during one month. As shown in Figure 11 the iso-variations between new and old bed levels, highlights a fairly large deposit in the dredged areas, estimated to have a maximum value of about 50 cm. This deposit contributes to sediment transport of particles from both banks of the estuary, causing a relatively important erosion which is estimated respectively at about 70, 50 and 40 cm at the level of the dredging site, the historical monument and the right bank of the beach. This erosion and deposition processes will continue in a non-linear tendency until an equilibrium state is reached where the local slopes become relatively low.

The compensatory solution
The fluvial sediment transport of particles which constituted the main sediment supply resource, are entirely stopped by the dam constructed in 1974. Therefore, the formation of a sand bar at the mouth of the Bouregreg estuary shelters the dredged areas and reduces sediment exchanges between the estuary and the marine  environment especially for fish activities. The supply of sediments to the dredged areas is therefore done, as confirmed by the mathematical model, by sediment transport of solid particles from the neighboring coasts. This erosion process is amplified by the combined action of the hydrodynamics of both tide and the dam's water releases and can endangers the environment and the stability of the neighboring monuments (Oudaya). To remedy to this critical situation, we recommend to control and well manage the dredging procedure in the estuary by avoiding intensive dredging in specific points.
For this purpose, the better solution choosed par decision-makers represent a dredging channel on the right bank in the northeast part of the es-tuary and passing through the bar ( Figure 12). This dredged channel, is extending from the old bridge to the piers at a depth of -3/ZH, contribute to improve navigation conditions in the estuary and to promote the exchange of water and sediments with the Atlantic Ocean. Also, subsequent dredging is needed to maintain sustainability conditions of the channel bathymetry.   Figures 13 and 14 show the simulations results. We can observe that in addition to the quantities of sand that will be dredged initially the maintenance of the channel will facilitate the circulation of water and especially guide the fluvial flow in the ebb phase, in order to compensate or reduce the negative impacts on the monuments. This is clearly confirmed by the model which shows a reduction of about 30% of erosion impact, observed before the digging of the channel (Figures 10, 11, 13 and 14). We can observe that this scenario makes possible the improvement of navigation conditions especially for medium and small boats, barks, water sports, etc. It permits also a better water circulation between the Bouregreg valley and Atlantic Ocean.

CONCLUSION
The Bouregreg valley is a complex system submitted to several negative impacts both under the anthropogenic and natural actions affecting continually water quality and morphodynamic processes. This paper assess hydrodynamic and sed-iment issues especially the sand barrier formation which is reducing exchange between the ocean and the estuary. A two horizontal numerical model HYSED is developed and implemented coupling hydrodynamic,sedimentary and morphological processes to study the impact of dredging on the estuary's morphodynamics. The model is used to assess the dredging operations impacts on the es-tuary with special emphasis on possible technical solutions to improve the situation. After calibrating the model by measurement data, multiple modeling. simulations for the dredged channel solution was realized. The results shows that in addition to the quantities of sand that will be dredged initially, the maintenance of the channel will facilitate the circulation of water and especially guide the fluvial flow in the ebb phase. The HYSED model results show a reduction of about 30% of erosion impact, observed before the digging of the channel. Therefore, this scenario makes possible the improvement of navigation conditions especially for medium and small boats, barks, water sports, etc. Also, it permits a better water circulation between the Bouregreg valley and Atlantic Ocean. However, there is a luck of numerical modeling of Moroccan's estuaries in general and particularly for the Bouregreg area. This issue is very important for a country like Morocco which have a strategic location on Atlantic and Mediterranean coasts (more than 3500 km of length); and a surface water resources management based on approximately 200 big dams constructed on the important rivers (Oum rebia, Moulouya, Loukos, etc.), with important impacts on the morphodynamic processes of the estuary's ecosystems. Therefore, there is an increasing need to improve cooperation between water authorities, estuary's managers, universities, users, consultancy and engineering professionals, to more understand the different processes involved in these complex ecosystems and to study and predict the main aspects involved like: salt intrusion, climate change impacts etc. The main objective of this cooperation is to better understand the estuary's processes and improve the efficiency of their management for a sustainable development tacking into account the environmental, economic and social aspects.