Numerical Modelling of Free Surface Agitation in a Coastal Lagoon by Roadway Path Influence

. The development and construction of highway infrastructure are essential in developing countries, whereas its layout and construction sometimes interact with the coastal environment. One problem to attend to is that the outline and geometry designs impact as little as possible on the hydrodynamic circulation of coastal bodies in order to alter the associated ecosystem as little as possible. The study area is located in the north of Colombia, and is made up of a continental coastal zone (Mallorquín Lagoon) and a marine zone (Caribbean Sea), in which a highway is projected that provides communication between two locations. This study presents the application of a numerical model previously developed and modi ﬁ ed by the Berkho ﬀ equation, which is developed in a ﬁ nite di ﬀ erence scheme and has been validated and applied in di ﬀ erent works in coastal and ﬂ uvial shallow water areas. The application of the model was carried out in a hydrodynamic circulation research project for a one-way highway through a coastal lagoon, where the knowledge of the magnitude of the incident wave height in the structure of the road body is necessary for the design, protection elements, and road geometry. Two numerical simulation scenarios were carried out, specifying normal conditions and extraordinary wave conditions in the month of November with a simulation time of 15 days, obtaining the velocity ﬁ eld associated with coastal currents, waves, and wave modi ﬁ cation phenomena, such as refraction, di ﬀ raction, and re ﬂ ection, which provide the height of the incident wave on the highway and the recirculation patterns in the coastal lagoon to identify alterations in the ecosystem. The results of the wave height in each scenario and the velocity ﬁ eld provide values to be used in the design, type of armor, and dimensions of the protection works required for the proper functioning of the road structure.


Introduction
The coastal physical processes consist of many systems and variables from different origins, such as the continental sources (rivers, estuaries, lagoons, and wetlands); the oceanographic ones (ocean currents, tides, waves, and coastal sediment transport); and the hydrometeorological ones (winds, precipitation, and temperature).The study of these ecosystems includes physical variables, which belong to different special and temporal scales.
Coastal aquatic ecosystems have some physical characteristics that have been studied by numerous scientists and whose results and conclusions have concluded that these natural systems are the most productive ones world-wide [1].The coastal lagoons are aquatic bodies, which most of them have permanent communication with the sea and they are the result of collisions of continental and saline water masses; likewise, coastal systems play a leading role in the life cycles of many marine species; on the other hand, the fishing activity in these areas is facing serious difficulties due to the lack of infrastructure for developing this activity and the lack of ecological regulation.
The structural element design in port works, such as dock, seawall, and breakwater, is mandatory to reinforce a region infrastructure; the sea behaviour, in coast zones, needs to be known and, hence, sea surge study becomes important.The processes of sea wave transformation include height variations and propagation direction, which are called refraction, diffraction, reflection, and breaking waves.
One of the tools that are contributing to previous hydrodynamic studies are computational numerical models applying fluid dynamics (CFD), which constitutes an approximation in the study of fluid dynamics equations that are used to describe different physical phenomena related to the movement of fluids, such as those produced by compressible and incompressible flows.The use of CFD allows different numerical calculations to be made under various conditions that form scenarios or case studies.
Specially, an analysis by numeric models is accomplished in this research, where-by, the behaviour of different hydraulic phenomenon can be predicted, they are important to estimate the projection and the planning of required port works.
The mild-slope equation (MSE) was developed since 1960 and currently modifications until today.In 1990s some software programs were created and commercialized to solve the numerical simulations in wave process analysis, especially studies in nearshore and wave interactions with coastal structures.The main goal was the presented a numerical model developed to solve free surface nearshore wave processes [2,3] and the interaction with breakwaters in construction on roadway.
The MSE describes the wave phenomena surface established by Berkhoff [4].The formulation analyses the refraction and diffraction effects mainly, whereas the wave breaking, bottom shear stress, and wave currents are calculated using a hydrodynamic model.
The numerical method used to solve the MSE was the finite-difference method [5], using alternating-direction implicit (ADI) method schemes [6,7].Its solution consisted in defined boundary conditions along the dry cells like marine structures and the requirement of solving a large number of simultaneous linear system of equations.Recently developments have been used bi-conjugate gradient method and implicit approximations [8].The use to finite difference methods on hyperbolic formulation were introduced by Song et al. [9] and Zhang et al. [10] with proposed to accelerate the computation and considerate three-point finite difference location.A light explicit finite difference method was presented by Zhang et al. [10] and otherwise an explicit scheme method by Lin [11] and Li and Fleming [12].
Actually, some research based on finite volume were presented by Inan and Balas [13], using curvilinear coordinates with ADI solution method [14] and spectral energy methods can be found at Tong et al. [15].Other research is available in Zhang and Li [16] with numerical solution of sea wave phenomena.
Some numerical computation limitations in elliptical approach MSE, is solve this formulation over a larger zone.To attend this problem, is used a set of three equations to perform the calculations of the phase velocity and wave height.
The numerical solution applied was a tridiagonal system to phase equations and pentadiagonal to wave height equation, this novelty in the implementation of a numerical solution allows solving the Berkhoff equation for waves in a simple way, optimizing computational resources, and giving an approximation of adequate solution values compared with other equations, such as Boussinesq used in commercial software packages, such as MIKE 21 and SMS, as well as free software packages, such as COULWAVE and FUNWAVE.
The numerical grid used is the staggered cell and developed coupled with a hydrodynamic numerical model [17][18][19][20][21] to study field velocity, bottom shear stress, and sediment transport in study area.The numerical model was proved with literature examples with analytical solution and de-fined boundary conditions, the results obtained are acceptable values compared with expected solution.
Since the key objective of this study is to determine the surge incidence magnitude and frequency in the interest area, it is important to consider different factors that produce the sea surface agitation.The principal agent responsible for generating agitation is the sea wave, which are produced by the wind friction on the aquatic surface.
The main oceanographic and meteorological data (wave direction and height, wind direction, and speed) were given by the Center for Oceanographic and Hydro-graphic Research of Colombia (CIOH by Spanish abbreviation) [22] and the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM by Spanish abbreviation).
Research consulting and analysis of key findings were necessary for the realization of this project.

Sea Surge Model Description. The equation developed by
Berkhoff was used, which is known as the soft slope equation.In addition, it is one of the suitable equations to simulate refraction-diffraction phenomena, in places where energy concentration, during the advance of the wave front, is produced by the seabed irregularities.Traditionally this equation is represented according to equation (1): Where the temporal change of the spectrum rate is represented by the first member of the equation, the energy propagation is contained by the second member, the wind entries by the third one, the sea-wave energy distribution among the nonlinear components by the fourth one, the 2 Modelling and Simulation in Engineering dissipation due to surge by fifth one, the losses caused by friction by sixth one and, finally, losses due to filtration are included by seventh one.
The numerical modelling is used to solve Equation (1), in an elliptic approximation (Equation ( 2)), which describes the propagation of a free and periodic wave on the surface and its finite amplitude based on complex bathymetry, where refraction-diffraction-reflection phenomena are represented by the deformation for wave approximation to shallow water and obstacles.
where η(x,y) is the level of the free-surface elevation (m), C(x, y) is the phase celerity or velocity (m/s), and C g (x, y) is the group celerity (m/s).The solution to the Berkhoff's equation is derived and expressed by the flow rate of the surge components X and Y [23].This was carried out by means of an implicit scheme in time by the following equations in finite differentiations.
where Q(x, y) is the flowrate in the horizontal plane and is given in m/s.
For the solution to the free-surface variation equation, due to the sea surge in function of previously mentioned flowrate, it is shown in the following equations: where h(x, y) is the group factor, k(x, y) is the number of wave (2π/L), and h(x, y) depth (m).This type of equation has the advantage of having the flowrate values, it and al-lows to find the sea surge direction and plan the presence of obstacles in a simple way.
The discretization of Equation ( 2) proceeds as follows: The solution of Equations ( 6) and ( 7) produces a tridiagonal system for every row i and j column, this implies that we must solve n i × n j lineal systems, for this reason it is important to have an efficient method, in this research we use the Thomas algorithm.6 Modelling and Simulation in Engineering The approximation for the wave height equation is given by: Equation ( 9) is a fourth order-accurate and produces a pentadiagonal banded lineal system.Such system requires the previous calculated phases (QX, QY), and the wave height in the previous time.We utilize a sparse matrix iterative method to solve the system.The values of η, QX, and QY, could be zero at time "t" except at the boundaries where the waves are incidents, this condition is necessary to determine the new values at time (t + Δt) [26].
For the wave height of the quasi-oscillating waves produced by the shock of the incident waves of height (H i ) versus the reflecting waves of height (H r ) the new wave height could be estimated as: The x-direction phases, the proper condition for Equation ( 7) at i = 1 and i = n i along the line (y = j × Δy + 1/2) is determined as: The results obtained are presented by Sato et al. [27] about wave height in the study of a physical model of a port (Figure 1), where the geometry of the port presents conditions of complicated borders that more closely resemble the reality of a port.
The incident period reported in the experiments was 0.73 seconds at an angle of incidence of wave to 0.0 °with respect to the axis Y m of Figure 1, besides considering that the breakwaters have a reflection coefficient of 0.8 that represents almost a null passage of the waves on the breakwater and a coefficient of 1 for open borders (total pass).
For the numerical model, the mesh is presented in Table 1 and the obtained results are presented in Figure 2.
A comparison of the qualitative type of the results obtained from the modelling is presented in Figure 3 with the Figure 1; analysing the resulting values of the wavelength model, they present a 5% decrease compared with those originally presented by Sato et al. [27], the reason was that a coefficient of 1 was used because the heights with the suggested coefficient they gave very small values to those presented by this author, later other reflection coefficients other than those originally used were tested to check the  Modelling and Simulation in Engineering behaviour of the model under different conditions reflection getting values totally different from the original ones, to which we conclude which is by the use of different computers.Finally, the computation time necessary to obtain the solution was 121 seconds to obtain a good approach.

Results and Discussion
The road structure project crosses in its route through a coastal lagoon known as Cienega Mallorquín, is shown in Figure 4 on the right margin, where it is located at the north of Colombia next to Barranquilla city and ends at the Caribbean Sea where a port zone will be built.
The main research goal is to determine the sea surge incidence and magnitude on the coastal marine and continental zone.Different factors were considered to stablish the generation of the sea surface agitation.Surge is the main agent, responsible for generating agitation that is produced by the wind friction on aquatic surface.
The roadway path along the coastal route is shown in Figure 4 on the right margin; in the same form, five control points are placed to monitor the wave height and estimate the level friction surface and roadway protection works (seawall) must be built.
The wave direction and height data were obtained from the CIOH.Ordinary and extraordinary sea surge conditions that were considered in the simulation.
Ordinary sea wave conditions to be considered in simulation: Direction: North.Wave height: 1.2 m.Period: 6.0 seconds.Simulation time: 15 days.Extraordinary sea wave conditions to be considered in simulation: Direction: North.Wave height: 3.6 m.Period: 8.1 seconds.Simulation time: 15 days.The tidal condition for both scenarios is shown in Figure 5, where the registration time is the complete month of November 2015.
The study domain was discretized in a finite difference grid (Figure 6), with spacing of 15 m × 15 m in each calculation cell.

Simulation without Landward Roadway
3.1.1.Normal or Ordinary Conditions.The results of the simulation analysis for November 2015 are presented in Figure 7.
The control points in the studied domain are shown in Figure 4b, which are located next to the roadway path.
Wave height results are registered in the control points, and they can be observed from Figure 8, where the sea free-surface evolution in time is seen in a 12,000-second simulation (less than one hour), and the wave height starts to get stabilized with surge periodic conditions.
3.1.2.Extraordinary Conditions.Simulation results for these conditions are not different from normal surge conditions, the main difference is found on the height that can be reached by the wave at a special condition.The sea-wave simulation in November is shown in Figure 9.
In the same situation, the wave heights registered at the control points are presented from Figure 10 for the extraordinary conditions.
A comparison with both stages from the three stations (1, 3, and 5) are presented in Tables 2 and 3, where the  3.2.Simulation with Landward Roadway.Consequently, simulation results for the area of interest with roadway included are now presented.For simulation, the roadway path structure was considered as a slope with a prefabricated concrete armor (dole or tetrapod) with a reflection coefficient of 0.12, so that the damping energy is considered at about 75%; the same simulation conditions are repeated with the same conditions and control points.Values registered at the control points are shown in Figure 12.The values registered in the control points with these conditions are shown in Figure 14.
Finally, a comparison between control points 1, 3, and 5 for both stages is presented in Tables 4 and 5, observing the changes among them not very noticeable, taking into account the surge elevation and the wave height for the month in simulation.
The results shown in the study allow us to observe the waves that appear on stage without the outline of the road in an open area, such as within the swamp; obtaining wave heights for normal conditions and extraordinary conditions (storm season), where the results of the numerical model are compared with measurements from the institutions in charge of measuring oceanographic variables, showing an acceptable approximation to the records and observations in place.The equations of the numerical model, previously developed by Herrera et al. [2], applied to infrastructure works, such as breakwaters, show an acceptable approximation to agitation phenomena in shallow areas, this was observed in the values obtained from the numerical simulation within the swamp and its interaction with the coastal zone.
The structure of the road projected in the marine area and within the swamp, is not yet constituted or defined in its elements or geometry, in the simulation an element with a permeability of 0.5 and a refractive coefficient of 0.7 was proposed, taking other structures as examples, existing coastal fences to the study area.
The results and precision of the wave height values, allow to have a preliminary design of the location of the line and possible geometry to reduce the effects of the incidence of the waves in the structure; likewise, mitigate secondary effects due to agitation within the swamp and, thereby, interfere with aspects of the marine environment.The numerical model presents acceptable results in a simulation or execution time, compared other complete and complex commercial models.

Conclusions
The wave heights in surge normal conditions with mean low-tide level in no more than 15 cm and a maximum of 0.44 m can be reached so the protection works and it must be designed for these sea wave heights.There is not any meaningful presence of agitation in relation with the wave height inside the coastal lagoon; for that reason, the protection works that cover the embankment road can be built of rock elements with average dimensions and weight, whose calculus can be accomplished by the Hudson's method (1956), [25] method.
A longer time study is recommended for simulation (about 6 months) to observe other meteorological phenomenon influences, such as hurricanes and rains that have big-ger impact on a region and increase the sea level and height and wave period over the influence of these values on the hydrodynamic field and its effects on the road structure; likewise, use a simulation time greater than 15 days.It is also suitable to observe the sediment transport to analyse the limitation line in the Cienega of Mallorquin, as well as the proposed roadway structure and its alteration in local coastal morphology.

Figure 1 :
Figure 1: Geometry and contour of wave heights within the port, Sato et al. [27].

Figure 8 :
Figure 8: Fluctuation of the free-surface at the control point 1, 3, and 5 in roadless simulation normal condition.

Figure 9 :
Figure 9: Development of wave propagation in the studied domain in extraordinary conditions: (a) 0.4 day (9.6 hour) simulation.(b) One day simulation.(c) Fifteen days simulation.

Figure 10 :
Figure 10: Fluctuation of the free-surface at the control point 1, 3, and 5 in roadless simulation extraordinary condition.

3. 2 . 1 .
Normal or Ordinary Conditions.Results are presented in Figure11, the same previous conditions for simulation were considered, with the only difference at observing the wave height behaviour by the roadway presence; a convenient stage for November is shown in the next figures.

Figure 11 :Figure 12 :
Figure 11: Development of wave propagation in the studied domain in normal conditions with landward roadway: (a) 0.4 day (9.6 hour) simulation.(b) One day simulation.(c) Fifteen days simulation.

Figure 13 :
Figure 13: Development of wave propagation in the studied domain in extraordinary conditions with landward roadway: (a) 0.4 day (9.6 hour) simulation.(b) One day simulation.(c) Fifteen days simulation.

Table 1 :
Parameters for the simulation of the incident wave in the study domain.

Table 2 :
Maximum wave height referred to below mean low water level for control point 1, 3, and 5.

Table 3 :
Average wave height referred to below mean low water level for viewers 1, 3, and 5.