A numerical model based on the mild-slope equation of water wave propagation over complicated bathymetry, taking into account the combined effects of refraction, diffraction, and reflection due to breakwater, is presented. The numerical method was developed using a split proposed version of the mild-slope equation and solved by an implicit method in a finite volume grid; this technique easily allows model the wave effects caused by the breakwater building in coastal waters, where industrial and other economic activities take place. Controlled case studies have been made and the results match very well with the reference solution. The capability and utility of the model for real coastal areas are illustrated by application to the breakwater of the Laguna Verde Nuclear Power Plant (LVNPP).

1. Introduction

Recently a great deal of public concern and interest are being directed towards the coastal waters, where industrial activities take place. As a result engineers are faced with and increasing demand for solutions to coastal problems. Among these is the understanding of the water wave transformation over the complicated bathymetry, from deep water to shallow sloping beaches, which dominates various aspects of coastal hydrodinamic processes, especially the generation of alongshore currents as well as sediment transport.

The mild-slope equation- (MSE-) based wave models emerged in late 1960s and have since reached total maturity by early 1990. A number of free and commercially MSE models are available today like CGWAVE and CREDIZ models, but in the field it is not possible to use them because the engineers do not have the enough time to create the meshes and wait for the results in field. The key of the model presented is the previous calculus of some important nearshore wave processes like wave current and wave-wave interactions made with a research hydrodinamic model early developed [1].

The mild-slope equation which describes the propagation of periodic, small amplitude surface gravity waves over a seabed of mild slope was first established by Berkhoff [2]. This original equation takes into account the combined effects of refraction, diffraction and reflection, while the influence of wave breaking, bottom friction, currents and wind is neglected; for this reason we utilize a hydrodinamic model [1] to calculate these effects.

The solution of the mild-slope equation was obtained by using the finite-element method [3, 4] and finite-difference schemes [5, 6]. However, this mild-slope equation is essentially of elliptic type and therefore its solution is very difficult due to the necessity of giving boundary conditions along the closed curve and the requirement of solving a large number of simultaneous equations. The latter restriction makes it computationally costly to find the solution for short waves over a larger area. To overcome this problem, in this research the mild-slope equation is split in tree equations to perform the calculations of the phase velocity and wave height.

The proposed numerical scheme requires the solution of a tridiagonal and pentadiagonal systems to solve the phases and the wave height, respectively. The type of grid used is the staggered cell with the possibility of a grid refinement in zones of interest. The numerical code developed was coupled with a hydrodynamic code [7] which provides the hydrodinamic in the study area. The numerical code was proved with examples where the analytical solution exists and has constant conditions obtaining acceptable values with an error less than 1%.

The study area corresponds to the breakwater of the Laguna Verde Nuclear Power Station located on the coast of the Gulf of Mexico, in Alto Lucero, Veracruz (LVNPP). The real meteorological parameters and the waves direction information of the area were acquired in order to obtain more accurate results.

2. Mathematical Model

The propagation of periodic, small amplitude, linear surface waves over a seabed of mild slope can be described by the following equation [8]: where is the wave height, is the phase velocity, and is the group velocity.

Equation (1) accounts for the combined effects of refraction, diffraction, and reflection and can be used for the computation of wave propagation in coastal regions and for the computation of wave transformation with a reflecting boundary such as harbour or breakwater [9]. However, it is essentially for elliptic form, and its solution needs a boundary condition along the closed curve and a great amount of computing time and storage. Therefore for short waves over an open coast area, it is infeasible.

The aim of the numerical solution proposed for (1) is to split the equation, in its phase velocities components and and the wave height . The equations for the phase velocities and in the directions and , respectively, are expressed as

The wave height is determined by where (group factor) is defined as follows: where denotes the wave number and the depth.

The advantage of dividing the general (1) in (2), (3), and (4) is that it models the presence of breakwaters in a simple manner.

3. Numerical Scheme

The numerical solution of (2), (3), and (4) is performed using a staggered grid illustrated in Figure 1, in the grid wave height is calculated at the solid grid points, labeled , , , , and , and so forth, and the velocities are calculated at the open grid points, labeled , , , and . Specifically, is calculated at points and and is calculated at different points and . The key feature here is that wave height is calculated at different grid points. In Figure 1, the open grid points are shown equidistant between the solid grid points, but this is not a necessity.

The discretization of (2) proceeds as follows: analogous for (3) the discretization is expressed as

The solution of (6) and (7) produces a tridiagonal system for every row and column, this implies that we have to solve lineal systems, for this reason it is important to have an efficient method, and in this research we use the Thomas algorithm. Figure 2 shows the computational molecule for the calculation of and .

The approximation for the wave height is given by The last (8) is fourth-order-accurate and produces a pentadiagonal banded lineal system. Such system requires the previous calculated phases ) and the wave height in the previous time. We utilize a sparse matrix iterative method to solve the system. The values of , , and could be zero at time except at the boundaries where the waves are incidents; this condition is necessary to determine the new values at time . Figure 3 shows the computational molecule for the calculation of the wave height.

4. Boundary Conditions

To obtain the wave propagation at open boundaries, we utilize Snell's law, which permits that the waves leave the domain without modifying their direction nor magnitude. The obstacles and buildings inside the domain are considered wall boundaries with a reflection factor , ranging from 0.0 (no reflection) to 1.0 (completely reflected).

To calculate the wave height at the borders of any breakwater, we formulate two expressions:

Equation (9) is used at the start of the breakwater of any length and (10) is used at the end. Figure 4 shows how the boundary conditions are applied.

4.1. Energy Dissipation

For the wave height of the quasioscillating waves produced by the shock of the incident waves of height versus the reflecting waves of height , the new wave height could be estimated as where is the factor of group in , factor of group in , is the distance with respect to -axis, and the distance with respect to -axis.

At the wall, the proper boundary condition for (6) at and along the line is determined as

Analogous for (7) at and along the line are used the next equations:

Finally, the angle frequency could be estimated by the relation between the phases and as

5. Controlled Numerical Experiments

To carry out the validation of the numerical code, theoretical examples were used in order to compare the results produced by the effects of refraction, diffraction, and reflection. In our first example we want to reproduce the diffraction and reflection caused by the presence of two symmetrical breakwaters with an opening in the middle.

Figure 5 depicts the geometry of the obstacles and the expected waves diagram at the exit of the opening.

In Figure 6 we can see the effects of diffraction and reflection, where the behaviour of the waves agrees with the expected solution depicted in Figure 5. The size of each grid cell is  m and the time step is of 1  seconds, the grid size is of and the breakwaters are located in the middle of the grid with an opening of 10 m.

In our next example we model the wave transformation produced by the breakwaters of a small rectangular harbour, depicted in Figure 7; the experimental solution (reference solution) for this problem was presented by Unluate and Mei [10] and a numerical solution by Lee [11] and J. Maa et al. [12].

The miniexperimental harbour geometry is of 0.3212 m in length × 0.0605 m of wide and 0.2576 m of depth, the angle of incident of the waves is zero degrees, the initial wave height is 0.01 m, and the initial parameters summary of the numerical test is depicted in Table 1.

Figure 8 shows the contour of the wave height obtained with the numerical code developed when the waves are stable; the boundary condition at the wall was considered with a reflection coefficient of 1.0 that means a total reflecting wall.

Later a quantitative comparison between the solutions obtained by Unluate and Mei [10] and J. Maa et al. [12] was made. Figure 9 shows normalized wave height obtained at the receiver compared with the reference solution and the numerical solution obtained by J. Maa et al. with a max relative error of 0.75%; the technique used to approximate the error could be seen in [13].

6. Numerical Experiments in the Breakwater of the LVNPP

The numerical simulation was carried out in the coastal area of the LVNPP which is located on the coast of the Gulf of Mexico, in Alto Lucero, Veracruz. It is the only electric power generating nuclear plant in Mexico and produces about 4.5% of the country's electrical energy. LVNPP has an original installed capacity of 1,365 Megawatts (Mw). It consists of two units General Electric Boiling Water Reactors (BWR-5) using Uranium (U235 Isotope 3% enriched) as fuel. Unit-1 (U-1) started its operation on July 29, 1990. Unit-2 (U-2) started its operation on April 10, 1995. The plant is owned and operated by the national electric company owned by the Mexican government (see Figure 10).

The generation of electric power at the plant is based on the technology of nuclear fission of uranium atoms, which takes place in the reactor. The energy released by the nuclear fission is transferred as heat from the fuel to the cooling water, which boils into steam. The water discharge of the cooling system is carried out in a trough channel with a flow rate of 63 m3/s and a mean speed of 1.4 m/s.

The type of the cooling system is considered open, actually the intake of the water is superficial, the breakwater which protects the intake has two piers, the main problem is the blockage, and the recirculation of hot water has causing deficiencies in the cooling system.

The grid of the study domain is shown in Figure 11, it is conformed by 5,896 cells with a constant size of  m.

The initial parameters for this simulation are defined in Table 2.

In Figure 12 is shown the result for the coastal area, we can see the refraction effects when the waves are near the coastal line produced by the bottom topography.

Studying in detail the behaviour of waves and their circulation inside the breakwater is very important to study the sedimentation effects, because the sediments transports have the potential to increase the maritime life producing the growing of different animals like mollusc, shells, and clams. This issue could cause the siltation and obstruction of the intake area reducing the performance of the cooling system.

Due to the importance of the breakwater area it is necessary to make a refinement of the grid into the breakwater; this refinement is made with 600 cells with a size of  m. In Figure 13 is shown the grid and the bathymetry of the breakwater area.

To probe the model, as first step, is proposed a breakwater geometrically idealized (see Figure 14); this means a controlled scenario to represent the wave transformations into the breakwater. This idealized breakwater is very close to the way and localization of the real breakwater to observe the behavior of wave propagation.

The results of the wave propagation are shown in Figure 15 for different simulation times; we can see that the diffraction and reflection phenomena are diminished by the breakwater structure.

Once the simulations over the idealized breakwater were done and show results with well-defined, patterns we carry out numerical experiment with the real bathymetry and the real breakwater form. The simulation result of the stable waves is depicted in Figure 16.

In Figure 17 we can see that the waves height into the breakwater is minimal; we find just a small disturbance in the free surface ranging from 0.05 to 0.1 m; this pattern is consistent with the average data acquired in the meteorological stations.

7. Final Remarks

In this research we develop a very fast numerical code for the study of the wave propagation over irregular bottom topographies. The model has produced consistent results with controlled cases and the result for the coastal area of the LVNPP is very acceptable compared with the data acquired in situ.

The aim of the model is to split the governing equation, in its phase velocities components and wave height, reducing significantly the computing time and the complexity of the implementation using sparse matrices. We can solve the system over a grid of 5896 cells in just 120 s using sparse matrix algorithms in a standard workstation.

The model is capable of reproducing the effects of the reflection and diffraction with accuracy and the refraction could be reproduced only for domains with cells greater or equal to 102 m. In future work we will present the sediment transport and particle transport coupled with the wave transformations.


:Phase velocity (m/s)
:Phase velocity in the direction (m/s)
:Phase velocity in the direction (m/s)
:Wave heights (m)
:Angle of frequency (angle in radians)
:Wave potential function
:Group factor
:Reflection factor.


The authors would like to thank engineer Ivan Campos for the information he provided to calibrate the numerical code. The authors would like to thank the reviewers for their helpful comments to improve the manuscript.