Three Dimensional Hydrodynamic Model with Multiquadtree Meshes

This study presents a three dimensional model for the transport of conservative contaminants, which can be used for bodies of water which are affected by winds and/or tides. The model solves the equation of mass transport, based on results obtained using a hydrodynamic model for shallow waters that works in a finite volume scheme and a type of hierarchical grid, called multiquadtree, which is adaptable to the bathymetry. To solve the vertical coordinates, the coordinate z is transformed into a sigma (σ) coordinate, thus allowing the same number of layers in the vertical, regardless of depth. This hydrodynamic model is validated using two cases: a long wave propagated in a channel of variable width and bottom and wind action in a rectangular basin. Finally, the results obtained are presented for a hypothetical single port outfall in the bay of Campeche, México. The model developed here is both quick and easy to use and is efficient when compared with models presented by other authors since it uses adaptable grids which allow detailed solutions to be obtained for areas of interest such as coastlines and the area around an outfall.


INTRODUCTION
The run-time and solution accuracy of a computational fluid mechanics model depends, essentially, on both the numerical scheme used and the quality of the mesh over which the governing equations and boundary conditions are discretized.
In this article a multi-quadtree, finite volume approach for modelling coastal hydrodynamics and water quality is presented. Compared to the use of finite difference methods, hydrodynamic modelling using a finite volume method appears, at first sight, to be more computationally intensive. However, in the multiquadtree case, it is seen that modelling is not only faster but also gives a more accurate representation of the hydrodynamics than finite difference methods using regular grids.
The multi-quadtree has been developed from the quadtree approach. This uses an adaptive mesh generation based upon iterative subdivision of a square domain to produce a hierarchical solution grid for the model. The resolution of the grid is determined by the gradients of the properties being represented, such as velocity, concentration of polluting agents and temperature. Where large gradients exist, the mesh density is greatest, whereas in areas where the gradients of the properties are small, the grid is coarse. This minimises unnecessary computational effort and reduces the time needed to run the model.
The quadtree technique has been applied in recent years to the study of an increasing number of fluid mechanics problems such as bifurcation in channels, river flow, coastal hydrodynamics, transport of pollution and so on, enabling improvements in accuracy to be made in shorter timescales.
The multi-quadtree presented here has the advantage over the quadtree in that it permits a subdivision of the entire modelling domain. This enables rectangular and other shaped domains to be modelled using a number of square quadtree subdomains. This new technique also has a simplified and more efficient numerical scheme.
The main features of the multi-quadtree approach are as follows: The generation of the multi-quadtree mesh is very fast: this minimises the need for extensive computational resources. The mesh is dynamic: it can be updated during modelling to give an improved representation of the phenomenon under study.
This type of mesh can be applied to any study region.
This chapter is divided into four sections. In the first one, the quadtree and multi-quadtree mesh generation is presented and the latter demonstrated with a case study. The second section describes the numerical model. The last two sections describe the validation of the model and its application to a study of the Campeche coast, Mexico.

HIERARCHICAL MESHES OF MULTI-QUADTREE TYPE
The generation of a hierarchical multi-quadtree mesh begins with the division of the entire region of study, such as a river or a coastal lagoon, into a number of adjoining square sub-regions. For each of these subregions a quadtree mesh is then generated.
Each sub-region is scaled to a unit square which is then iteratively subdivided into 4 sub-squares. This continues at each locality in the domain until either the maximum specified level of mesh density is reached or else specific criteria are locally satisfied. These criteria are derived from quantities such as the local depth, vorticity or velocity fields. In this way the quadtree mesh is produced.
The generation of a sub-region quadtree mesh requires 2 steps: initial generation and regularization. Ancillary procedures to assign the geometrical and topological descriptions of the grid are also performed. After the generation of the sub-region quadtrees a third step, verification, is then necessary in order to ensure that the sub-regions are mutually compatible, leaving the final multi-quadtree mesh.
There are several examples of quadtree generators which may be found in Greaves [1] , Borthwick et al. [2] , Rogers [3] and Stallard [4] . The generator used in this article is based on that proposed by Stallard [4] . The improvements and additions are described in the following sections.
It might also be noted at this stage that whilst the quadtree technique is used to solve problems in 2D, it can be extended to 3D. In this case an octree mesh could be used, based on the iterative sub-division of a unit cube into 8 smaller cubes.

Initial generation:
The objective of this process is to produce a mesh of variable density by subdividing the region according to the initial data or seed points. These initial data comprise the following three types. For each of the m n sub-regions the seed points are plotted and the sub-region is scaled to a unit square. Then, the first two seed points are considered and the iterative division of the sub-region proceeds until either of the following criteria is satisfied: The maximum level of subdivision is reached. Only one single point exists inside any cell.
Ubsequent pairs of seed points are then taken and the process of subdivision continues for all the points in the sub-region. Figure 1 shows an example quadtree mesh generated around 4 points.
During the model runs, the level of subdivision may also be dictated by the magnitude of variables in that locality. That is, if during the modelling a parameter value exceeds a predetermined threshold, then the mesh density is automatically set to the maximum level of subdivision.
For each cell of the final mesh the following data are obtained: cell reference (identification number), cell coordinates, cell width, level of subdivision and also the coordinates and level of subdivision of the neighbouring cells.

Cell Numeration:
The subdivision of a cell consists of the creation of four new cells that represent the four To each daughter cell a positional index N m is assigned for each level of subdivision, m, according to the convention in Fig.2 The complete cell reference, N, for the cell in Fig. 1 that contains point A is thus 1411 and for point D, 1423. The advantage of this over other Quadtree generators is that it more compactly describes the cell location. That is only a single number is required for each daughter cell in the range 1-4 and not coordinate numbers. Posada et al. [5] . Regularization: The objective of regularization is to balance the model efficiency and mesh uniformity by prohibiting rapid spatial variations in cell size. In practice, the mesh obtained after initial generation is re-adjusted so that neighbouring cell dimensions do not differ by more than a factor of two, i.e. the difference in subdivision of neighbouring cells cannot be more than one level. In this way, the need for complex interpolation is avoided.
The process begins by selecting a cell that has the highest level of subdivision. The neighbouring cells are then adjusted, as required, until they are within one level of subdivision. The process is then continued for the neighbours of the neighbours of each cell.

Verification:
The third step for the creation of a multiquadtree mesh is the verification that all the sub-regions have the same cells in common. For this procedure, a new list is created of the coordinates for the centres of all the cells that form the boundaries of two subregions. These points are added to the original list of points, then the duplicated points in this combined list are erased and the whole process is repeated until the sub-regions agree.

Summary of advantages:
The advantages of the quadtree and multi-quadtree approaches are: Despite the apparent complexity of the resulting meshes, the generation of the mesh is automated. The structure of the data for the flow variables is arranged hierarchically. The ease with which the method can be refined and adapted for areas where the local gradients of quantities are strong. No transformation of coordinate systems is required The multi-quadtree can be used to model nonsquare domains, it provides improved efficiency for modelling large rectangular domains. However, the authors have found that if the relation between sides in a rectangular cell is greater than 3, numerical problems may be appear.
The main feature of the multi-quadtree method is that the mesh can be repeatedly and rapidly updated according to pre-established criteria derived from flow quantities such as vorticity, velocity and depth. The calculation time for refining and updating the mesh is insignificant compared with the overall modelling time, yet significantly improves the model performance.
Interpolation: Following the generation of a multiquadtree mesh many adjoining cells will be of different sizes. The regularisation process ensures that the difference in size is limited to a factor of two. This results in a virtual node on the boundary of the larger cell. The use of straightforward linear interpolation for the creation of this node has been shown to provide sufficient model accuracy Rogers [6] .

MULTI-QUADTREE MESH EXAMPLE: ACAPULCO BAY
By way of an example, the multi-quadtree mesh generation is now applied to the Bay of Acapulco, México, (Fig. 3). For this illustration the study region is first divided into 2 2 square sub-regions. These are numbered using the father-daughter convention (Fig. 2). Figure 4 shows how the process of verification regularizes the mesh across the boundaries of the four For example, note how the mesh at the centre of the entire region has been adjusted.

2D AND 3D NUMERICAL MODELS AND CONTAMINANT TRANSPORT MODEL
In this article, the hydrodynamic behavior of a body of water is obtained using the technique of splitting, which allows the division of a 3D model into two sub-models. In the first one, the depth averaged equations are resolved in the vertical, obtaining the mean velocities in the X and Y directions for each cell. The velocity in the X direction is known as U and that of the Y direction, V. The variation of the free surface in relation to the mean sea level is known as . In the second sub-model the results of the 2D model are considered as the initial data and from these the velocities for each cell in the x, y and z directions, referred to as u, v and w, are obtained.
Background: Amongst the numerical schemes that can be implemented to solve the non-linear shallow water equations there are three common approaches: finite differences; finite element; finite volume. In recent years, the finite volume method has been used more frequently due to improvements in computational resources and its superior handling of discontinuities and flows in closed boundaries.

2D model:
The numerical model solves the governing equations with a finite volume scheme in two dimensions. This is based on a Godunov scheme using the hierarchical mesh. To solve in-viscid flows we employ a Riemman solver to calculate the approximation proposed by Roe, in a similar way to Posada [7] . The time integration is carried out with a first-order Adams-Bashforth technique.
The depth-averaged equations, solved by the model in 2 dimensions are: Where, H = h + = Free surface elevation with respect to the mean water level. h = The bottom depth with respect to the mean water level.

Momentum equation in X-direction:
Where C d is the drag coefficient of the wind, which has a value of 0.0026 [8] . The shear stress due to friction at the sea bed is calculated with the expression where, C D is the friction coefficient, which is a function of the Chezy coefficient

Eddy viscosity coefficient in 2D:
The numerical model provides two options. The first is to consider this coefficient a constant. The second is to consider the viscosity coefficient as a variable, calculated with the following expression.
3D model: The 3D model works on sigma coordinates (σ) to represent the vertical coordinate. The horizontal gradients are calculated using Cartesian coordinates [10] . Working with sigma coordinates allows the model to have the same number of layers, regardless of water depth. This is not the case with the 3D models which work with Cartesian coordinates, so that in deep waters there are a greater number of layers than in shallow waters.
In order to use the sigma coordinates (σ) it is necessary to transform the equations that govern the flow. This procedure is given in detail in Posada et al. [7] . The sigma transformation is defined as: There z = depth of the layer with respect to the mean sea level. As can be seen, varies between 0.0 and -1.0, 0.0 for the free surface and -1.0 for the bottom, as observed in Eq. 11. The transformation of the equations is made with the following expressions: There: x x *; y y*; t t *; The continuity equation, transformed into sigma coordinates (σ) is: The equations of angular momentum in directions X and Y are respectively: x y x y y These three-dimensional equations are solved with an algorithm in finite differences which calculates, implicitly, the vertical diffusion with a Crank-Nicholson scheme. The vertical convective acceleration is calculated with a first order upwind scheme, as proposed by Castanedo [10] .
From the equations of angular momentum in directions X and Y the total provisional velocities, i ,k (1) u , are obtained. To assure the correct connection between the two and three dimensional models, the vertical variations of the velocity with respect to the depthaveraged value, u ' i,k , are calculated with the following expression: Finally definitive velocities are calculated with the next expression: Where, U i = Depth-averaged velocity calculated with the 2D model. Figure 5 shows the decomposition scheme of variables. Once the horizontal velocities and for each layer are obtained, the vertical velocity is calculated from the continuity equation.

Coefficient of eddy viscosity:
In this article a coefficient of eddy viscosity is adopted with a parabolic vertical profile. This coefficient is calculated with the following formula: ( ) There G(n,z,t) considers the discharge of a pollutant agent in cell n and layer z.

VALIDATION OF THE HYDRODYNAMIC MODEL
In this section the hydrodynamic numerical model is validated with analytical solutions and compared with other published results. The 2D model is compared with the analytical solution of a long wave propagating through a channel with variable width and depth. The 3D model is compared with the solution by Koçygit [11] , where the effect of the wind in a rectangular basin is evaluated.
Propagation of a long wave in a channel of variable width and depth: For this case a long wave with amplitude of 0.02 m and period of 100 sec, is generated in a channel with the following dimensions: length 2500 m, width b = nx, where n = 0.2 and the initial depth at the channel mouth is 25 m. The depth, h, varies along the length of the channel according to, h = mx, where m = 0.01. The analytical solution presented by Rahman [12] is used. This does not include dissipative terms so these are neglected.
Bathymetry: For this validation, the length of the channel is 4 times greater than its width. Thus a multiquadtree mesh of 4 by 1 was specified with the seed point depths in the range from 25.0 to -0.1 m. Figure 6  The generated mesh (Fig. 7) has the following characteristics: The four sub-meshes (left to right) have 982, 826, 634 and 430 cells respectively. The smallest cells are 19.6 19.6 m, whereas the largest are 76.8 76.8 m.
Results: Figure 8 shows the variation of the free surface of the water in the channel at maximum level. These extreme values are reached where the channel is narrow and shallow. The maximum and minimum free surface values obtained with the numerical model are 0.181 and -0.226 m. Note that 3 tidal waves fit in the channel.
The velocity field in each of the 4 quadtree submeshes are shown in Fig. 9. These results show the point of maximum free surface height, for which a corresponding maximum velocity of 0.30 m sec 1 is observed. At the point of minimum surface elevation the maximum velocity observed is 1.0 m sec 1 . It should be mentioned that the velocity is greater in the case of the minimum free surface elevation since the water depth of this cell is small. In fact this cell is almost dry.
In each sub-mesh the multi-quadtree method dictates that the flow speeds are normalized. Thus in Fig. 9 the velocity scale varies between sub-meshes (as indicated by the labelled vector above each plot). These different scales have been retained for the purpose of visualization.
Analytical solution: Rahman [12] proposes the following expression for calculating the free surface in a channel with variable width and depth Here J 1 is the first order Bessel function. Comparison with analytical solution: Figure 10 presents a comparison of numerical results with the analytical solution. In this figure the effect of three waves is observed. In particular, it can be seen that in the first half of the channel both solutions are practically identical. When comparing the numerical maximum free surface elevation with the analytical form, its amplitude is slightly smaller.

3D circulation induced by wind:
In order to illustrate the circulation induced by wind in a rectangular lake a comparison with the example proposed by Koçyigit [11] Figure 11 and 12 show the longitudinal pattern flows. It is observed that the surface velocities follow the wind direction while at the bottom they go in the opposite direction. At the borders the flow changes direction 90 degrees to follow, first, the contour of the wall and later the bottom. This pattern agrees with that obtained by Koçyigit [11] , but with a better representation of the flow close to the walls, in the model of Koçyigit [11] velocity vectors jump at bottom, whereas in this model the transition is gradual.

APPLICATION TO REAL BATHYMETRY
In this section the 3D hydrodynamic model is applied to a real case study: Campeche Bay, in the state of Campeche, Mexico. The objective of this case is to see whether the artificial port of El Embutido acts as a barrier between a contaminant discharge to the west and the protected area to the east. Figure 13 shows the bathymetry which covers an area 25.5 km in length and 12.5 km in width. The maximum depth of 5.0 m is located in the northwest corner. The coastline runs between the port known as El Embutido and the Mexican Naval base of Lerma, indicated by means of the blue ellipse in Fig. 13.

Multi-quadtree mesh:
The model domain was subdivided longitudinally into 2 meshes (Fig. 14). These are labelled, West to East, as sub-regions 1 and 2 with 2173 cells and 2389 cells, respectively. The maximum spatial resolution was required between depths of 0.80 m and -1. Hydrodynamic modelling: In order to study the hydrodynamic behavior of Campeche Bay an astronomical tide of amplitude 0.65 m incident on the northern boundary borders of each one of the meshes was modeled. The following results correspond to the fifth tide cycle obtained by the bidimensional model, when the solution has already stabilized.

High tide:
The maximum elevation appears in mesh 1 and is equivalent to 0.28 m. The highest value in mesh 2 is 0.2 m. Next, the velocity fields in the two sub-meshes are shown, (Fig. 15) Maximum bidimensional speeds are found in mesh 2, equivalent to  Figure 15 shows the corresponding velocity field at low tide, again for one fifth of the tidal cycle. The velocity vectors are opposite to those obtained at high tide. The maximum velocity of 0.4 m sec 1 occurs in sub-mesh 2.
The velocity values calculated with the numerical model are coherent with the process modelled. In high tide the vectors move inland, while in low tide the direction is seaward, with maximum values near to the coast, in the shallow waters of the domain.  Figure 16 presents the hydrodynamic behavior obtained by the threedimensional model for three different layers, (upper, middle and bottom) at high tide. The pattern of circulation in all the layers is similar to that obtained by the bidimensional model, but magnitudes of vectors are not. For upper layers, these are greater than the average speed, while the lower ones are of smaller magnitude, maintaining the same direction.
In layer 01 the maximum velocity calculated is close to 0.60 m sec 1 , in the intermediate layer, 0.40 m sec 1 , while for layer 5 the speed is 75% less than that calculated in layer 01, due to the effect of friction.
These results are shown in Fig. 17 in which the results obtained with the models in two and three dimensions are compared. The corresponding cell is 1500 m off the beach in mesh 2. The bidimensional model predicts an average velocity of 38 cm sec 1 . In the three-dimensional model for layer 01 the velocity is 57 cm sec 1 and for the layer nearest the bottom, 0.10 cm sec 1 .  Figure 18-20, show the behavior of the polluting agent (mg L 1 ) for the high tide of the sixth tidal cycle, 15 h after the beginning of the discharge. Analysing the distribution of the contaminant it is concluded that the port of El Embutido is not an efficient barrier; in all the layers significant concentrations of the polluting agent are present near the beach and on the rotected eastern side.

CONCLUSIONS
Mesh system, based on the quadtree technique has been presented. This allows a region to be subdivided into several square sub-meshes. This allows non-square domains to be modeled and also improves numerical efficiency, since calculation time increases non-linearly with the mesh dimensions, providing considerable savings in calculation time over regular meshes. Furthermore, as the mesh is quick to generate, adaptive meshes able to update between model steps can be employed at little time cost.
The numerical solutions obtained with quadtree meshes are of similar accuracy to regular meshes and enable the use of the most sophisticated finite volume modeling approach. However, one of the disadvantages of the approach presented here is that the solutions are not always well-matched at boundaries where the mesh is coarse.
To demonstrate the possibilities of working with multi-quadtree meshes, a 3D finite volume hydrodynamic model was developed from the shallow water equations.
This model was shown to be fast and stable.
The model was validated with two idealized cases: firstly, the case of a long wave propagation in a channel of variable width and depth; second, the case of wind driven currents in an enclosed rectangular basin. In both cases, numerical results were in good agreement with the analytical and other numerical solutions.
Finally the hydrodynamic behavior of Campeche Bay response to an astronomical tide was presented. The numerical model predicted amplitude for this tide of 0.65 m at the coast, which agrees well with the field measured value of 0.70 m.