Three-dimensional hydrodynamic model for wind-driven circulation

This paper presents the development of a three-dimensional hydrodynamic model for wind-driven circulation using the moving element method in the vertical discretization. As most three-dimensional models use sigma coordinate transformation in the vertical discretization, the use of the moving element method is an alternative to the latter, with the advantage that there is no need for fixed subdivision of the water column. In the proposed model, the coupling of two- and three-dimensional hydrodynamic model is considered. The shallow-water equations are integrated in the vertical direction, finite elements are employed in the spatial discretization, and finite differences in the time discretization, to resolve the position of the free surface (ζ), and the vertically integrated velocities components. The three-dimensional module is used to compute the velocity profiles. In the three-dimensional module is employed the moving element method in the vertical discretization. The efficiency of the model is demonstrated through the comparison of its results with laboratory experimental results and with another model that uses sigma coordinate transformation in the vertical discretization.


Introduction
Economic and social development in areas close to water bodies has caused growing environmental risk, not only in coastal regions, but over the entire supply system: this imbalance is caused fundamentally by the release of industrial and domestic waste and the disorganized natural resources exploitation in these areas. For this reason, in recent years, the environment preservation has become a relevant factor in promoting social wellbeing.
One of the main supporting tools in environmental management is the numerical models such as hydrodynamic models. The characterization of the hydrodynamic circulation of coastal water bodies by numerical modeling is essential to assess and to monitor environmental impacts suffered in these regions; in this way, numerical modeling is becoming an indispensable environmental management tool.
The three-dimensional mathematical models, describing the hydrodynamics circulation of the water bodies, are well known, and are formed by the continuity equation and the Navier-Stokes equations. The numerical methods generally used to solve this system of equations are the finitedifference (FDM), finite-element (FEM), and finite-volume methods (FVM). The FDM [2,10,12,25] has been widely used in the area of fluid flow, applied to uniform grids. A significant advantage of the FEM [14,16,26] is its easiness in describing complex geometries. Recently, the FVM [4,6,9,13] has been successfully applied in computational fluid mechanics because of the easiness of association between physical and mathematical models.
The multilayer model is a three-dimensional model wide-ranging in the literature [15,29], in which the water body is divided into homogeneous layers solved according to a two-dimensional model, with integrated equations Technical Editor: André Cavalieri.
& C. L. N. Cunha cynara@ufpr.br along the vertical direction. Kinematic and dynamic boundary conditions are imposed along the interfaces. Another usually used three-dimensional model is known as the Quasi-3D (Q3D) [10,11,13], in which the equations are simplified by adopting shallow-water conditions, that is, by applying hydrostatic distributions for pressure fields. In this model, vertically integrated shallow-water equations are used to solve the mean velocity field, and for every time step, the velocity profiles are determined at each node of the grid. Another three-dimensional models have been developed, that solve vertical and horizontal modules separately, e.g., [5,17,19,32]. In these models, the free surface position (elevation) and the components of vertically integrated velocities are determined using two-dimensional equations, forming the two-dimensional module (2DH). In the three-dimensional module, velocity fields are solved based on three-dimensional momentum conservation equations, since only the derivatives of each variable in the vertical direction are implicit. An additional model can also be included to determine the density field, which implies calculating the salinity and/or temperature transport equation.
Most the three-dimensional models use sigma coordinate transformation in the vertical discretization. The model presented in this study employs the moving element method in the vertical discretization [24] as an alternative to the sigma coordinate transformation, since there is no need for fixed subdivision of the water column.
The moving element method has been proposed in many different contexts over the last decades. However, in most of the applications, the moving element method was used only to solve the bi-dimensional flow. Scudelari [24] used the moving element method for the solution of the bidimensional shallow-water equations and compared its results against the analytical solution, using irregular mesh, obtaining a good agreement. Stockstill [27] describes a method for determining implicitly the waterline and flow variables in shallow water. The model uses an implicit Petrov-Galerkin moving-finite-element representation of the two-dimensional shallow-water equations. The model offers a viable means of representing shallow-water flows where the boundary locations are not known a priori. Baines et al. [1] used a moving mesh finite-element algorithm for the adaptive solution of non-linear diffusion equations with moving boundaries in one and two dimensions, applied in problems involving time-dependent partial differential equations with moving boundaries. Vanzo et al. [28] derived a novel numerical method to discretize the system of governing equations for the pollutant transport by shallow-water flows over non-flat topography and anisotropic diffusion. The resulting method is implemented on unstructured meshes, and is assessed for accuracy, robustness, and efficiency on test problems including non-flat topography. The novelty of this work is the use of the moving element method in the vertical discretization. The advantage of this method is the versatility of the discretization, once there is no need of a fixed water column division, which allows a better adjustment in the vertical discretization in the water bodies with non-flat topography.
The three-dimensional model, presented in this study, is calculated in two modules: in the 2DH module, the free surface position (f) and the components of vertically integrated velocities (U and V) are calculated using the continuity equation and momentum conservation equations in the x-and y-directions, vertically integrated. The decoupling method employed to solve the two-dimensional model prevents the simultaneous solution of several unknowns in the system of equations [7]. Thus, only the free surface position is implicitly solved, while velocity components are explicitly solved. The three-dimensional model uses three-dimensional momentum equations, in the x-and y-directions, and the three-dimensional continuity equation to calculate velocity profiles. This model applies the moving element method in the equations discretization and does not use sigma coordinate transformation. The technique allows greater versatility in the discretization of water bodies exhibiting large bathymetry variations, creating a vertical discretization that is updated at each time step. This discretization is based on the grid used in the 2DH module. The two modules are coupled using free surface gradients and bottom shear stress. The present study describes only the three-dimensional module. For additional information concerning the two-dimensional module, the reader is referred to Cunha and Rosman [7]. The model was tested against laboratory data for windinduced currents, developed by Yu [31], and compared with another three-dimensional model that uses sigma coordinate in the vertical discretization [23].
Rosso and Rosman [23] developed a three-dimensional model for the Navier-Stokes equations with shallow-water approximations, considering hydrostatic pressure approximation, named FIST3D. In the FIST3D, the vertical spatial discretization is carried out through finite differences using sigma coordinate transformation, whereas the spatial discretization in the horizontal x-y plane is done through subparametric Lagrangian finite elements. The model is composed of two modules: a depth-averaged or 2DH module, through which the free surface elevation and the vertically averaged 2DH current velocities are computed; and a 3D module, that computes the vertical profiles of the horizontal velocity field, as well as the vertical component of the velocity vectors. In FIST3D, the turbulent stresses and turbulent viscosity coefficient parameterizations are identical to the model presented in this paper.

Mathematical model
The three-dimensional model is calculated in two modules: first, the depth-averaged shallow-water equations are computed: the free surface position (f) and integrated velocity components in the vertical direction (U and V) are calculated [7]. Then, in the second module, the vertical distribution of velocity components is determined. The basic equations used in this module in Cartesian coordinate system (x, y, z), aligned to the east, north, and vertical directions, are as follows: Equation of continuity: Equations of motion: where the unknowns are u, v, and w, representing the velocity components in directions x, y, and z; g is gravity; q o is the specific reference mass; t is time; s ij (i = 1, 2, 3; j = 1, 2, 3) corresponds to turbulent stresses, and f is the Coriolis factor, which represents horizontal acceleration, denoted by:2Xsenh, where X is the angular velocity of the Earth and h is the latitude of the site considered. In this case, the effects of vertical Coriolis acceleration are disregarded. Turbulent stresses can be written in a parametric form by following the approach presented by Rosman and Gobbi [22], which is based on the filtering techniques [21]. These stresses are written as follows: where i, j = 1, 2, 3 and k = 1, 2, 3, 4 with k = 4 corresponding to time t (in this context, x 4 = t), and K ij are the turbulent viscosity coefficients. The spatial and temporal Gaussian filters widths, in the x k dimension, are defined as K k = a k D xk , where a k is a homogeneous scaling parameter in the x k dimension. The value of a k calibrates the filtering terms magnitude. Usual values of a k are between 0.25 and 2.0, and most often good results are obtained with a k = 1.0. It is easy to verify that the filtering terms, as displayed in Eq. (4), behave as self-adjusting sub-grid scale turbulent stresses. It is also verified that in a non-structured finite-element discretization mesh, the magnitude of these terms is a function of the local resolvable scale.

Boundary conditions
When the vertical distribution of velocity components is determined, boundary conditions must be specified at the free surface and the bottom.

Bottom boundary condition
For the boundary condition, bottom velocity is equal to zero, that is: Coupling between the three-dimensional and two-dimensional models is achieved through bottom shear stress. Thus, bottom shear stress is defined as follows: where u * is the characteristic friction velocity, obtained from the profile of velocities calculated in the 3D model, C is the Chèzy coefficient, and U i represents the integrated velocity component in x i -direction, previously calculated in the two-dimensional model (2DH).

Free surface boundary condition
The boundary condition is the specification of surface stresses, in directions x and y. These stresses can be parameterized as follows: where q ar is the specific mass of the air, W 10 represents the wind velocity, measured 10 m from the free surface, b is the angle of wind direction in relation to the x i -axis, and C D is the drag coefficient, defined by Wu [30] as follows:

Land boundary conditions
For land boundary points exhibiting significant inflows or outflows such as a river or a small estuary ending in a bay, for an instantaneous flow situation, prescribing only the normal velocity component is not sufficient to define a well-posed problem. In fact, when boundary segments have significant inflows, in addition to the normal flux or velocity component, the tangential flux or velocity along the segment must be specified, and it is set equal to zero.

Parameterization of the coefficient of turbulent viscosity
As filtration represents the greatest portion of horizontal momentum diffusion, the horizontal turbulent viscosity coefficient has a little influence. By contrast, vertical turbulent viscosity coefficients are more important and essential to the predictive of numerical models capacity. The vertical turbulent viscosity coefficients follow parabolic distribution and can be defined, according to Fischer et al. [8], as follows: is the total water depth, and a is a calibration parameter.
The characteristic friction velocity can be calculated based on stress at the surface (s S i ) and at the bottom (s F i ), or expressed as a function of mean flow variables. Considered as a function of surface and bottom shear stress, characteristic friction velocity can be defined by the expression below: where u S Ã i corresponds to the surface characteristic friction velocity and u F Ã i is the bottom characteristic friction velocity.
Bottom characteristic friction velocity is assumed to be a function of the logarithmic velocity profile. In other words, close to the wall, velocity distribution can be well represented by a logarithmic function as follows: where u i (Dz) depicts velocity at the first level above the bottom and e is the length of equivalent bottom roughness. The above formulations show that the characteristic friction velocity can take on several values, which may differ substantially amongst themselves. Jin [11] concluded that the best determination for u * follows the choice between maximum values of u S For flows where bathymetry changes abruptly, the nearbottom velocity profile is subject to the effects of recirculating currents. Selecting characteristic friction velocity in accordance with Eq. (12) has a little significance for flow, given that the velocity profile tends to invert close to the bottom. In such cases, another characteristic friction velocity can be calculated, which is based on mean flow variables, by the coupling of the 3D and 2DH modules. Characteristic friction velocity for 2DH is denoted by the following: The characteristic friction velocity is then determined according to the expression below: As previously mentioned, filtration represents a substantial portion of horizontal momentum diffusion. Therefore, a simplified formulation for the turbulent viscosity coefficients can be adopted, which is obtained from the vertical means of Eq. (10), considering j = 0.4. Thus, the horizontal turbulent viscosity coefficients for the two-dimensional model can be defined as follows: where K V ij are the diffusion coefficients for conservation of momentum due to vertical averaging and K H ij are the diffusion coefficients for horizontal momentum conservation. The model solves the vertical and horizontal modules separately; the vertical turbulent viscosity coefficients are depicted in Eq. (9) and the horizontal mixing coefficients are shown in Eq. (15).

Numerical model
The numerical model uses a combined system: the implicit factored scheme in the time discretization (finite differences), the finite elements in the two-dimensional spatial discretization, and the moving element in the three-dimensional module spatial discretization. In both the twoand three-dimensional modules, the numerical method uses a decoupling technique, where only one variable is calculated implicitly and the remaining variables are determined explicitly. This decoupling is carried out via the Successive Substitution Method (SSM).
The calculation process employed for the three-dimensional model, denominated MOTRID, can be described as follows: the free surface position (f) and the vertically integrated velocities (U and V) are calculated in the 2DH module.
With the free surface position calculated for the entire domain, the three-dimensional velocity components u and v are obtained from the equations of motion (Eqs. 2 and 3) and w, from the continuity equation (Eq. 1).
Since the free surface position is calculated using vertically integrated velocities, three-dimensional velocity profiles must be corrected to satisfy the equation of continuity, such that the vertically integrated velocities in both modules are equal. With this purpose, initially, the values obtained from the three-dimensional profiles integration are compared with those arising from the two-dimensional module; afterwards, the velocity profile is adjusted so as to produce the same vertically integrated velocity of the twodimensional module.

Temporal discretization via implicit factoring
The reference system used in the three-dimensional module is shown in Fig. 1, where the origin of the z-axis is considered on the free surface.
In order to uncouple the equations, the velocity components u and v are explicitly written in the momentum equations as functions of the free surface position (f) by making use of extrapolations and interpolations. The following notation is employed in the discretized version of Eqs. For the extrapolated or interpolated variables, with a second-order scheme, the following quadratic approximations with three time levels are used [20]: Using the approximations defined in the implicit factoring scheme, the time discrete equation of motion in direction x is given by the following: where: The values of u 1 in Eq. (18) are as follows: and Eq. (19) can be rewritten as follows: where Similarly, the value of v 1 is calculated according to the following equation: where: The value of w 1 is determined from the continuity equation. The u 1 and v 1 values employed to calculate w 1 are not the extrapolated values, but rather those corrected according to the velocity profile, as mentioned before. This calculation is performed as follows: The boundary condition for Eq. (22) is as follows [18]:

Spatial discretization
The moving element method, which combines finite-difference discretization with interpolation functions typical of the finite-element method, was adopted for the spatial discretization of the three-dimensional model. The elements, for which interpolation functions are constructed, are defined for each discrete point on the grid, see Fig. 2, and their connectivities are based on the two-dimensional grid, see Table 1. The selected element is an 8-node hexahedron: both the unknowns and geometry are approximated by linear interpolation functions, see Fig. 3. The set of linear interpolation functions for the geometry of the element and the remaining unknowns are denoted as follows [3]: Considering the framework shown in Fig. 2, the twodimensional grid defined on the plane of vectors i and j, where the values of the free surface position f ? (x,y,t) are calculated, produces a three-dimensional grid where, for each vertical, the equations of motion and the continuity equation are solved separately.
The derivatives involving horizontal gradients can be defined for any point of the element provided that the global coordinates (x o , y o , z o ) and, consequently, the natural coordinates (n o , g o , s o ), are known. Using a generic variable f [x o (n, g, s), y o (n, g, s), z o (n, g, s)], its derivatives can be calculated as follows: where J T is as follows: Vertical derivatives are approximated using the finitedifference scheme. The framework depicted in Fig. 4 is used to represent the vertical system.
With respect to Eq. (20), Q x,k and P x,k are known coefficients, obtained from extrapolated/interpolated velocities and values previously calculated in the twodimensional module, such as the free surface position, f ? (x, y, t). By assuming ow oz % 0 and K xz ¼ m x , one has: where: The vertical derivates are approximated according to the following: Equation (28) can be written as follows: in which: The resulting system of equations, which already includes bottom and surface boundary conditions, is solved using the TriDiagonal Matrix Algorithm (TDMA). Note that since this model was developed for tridiagonal matrices, some values must be extrapolated.
Similarly, the velocity component in y-direction is calculated using the momentum conservation equation at y, while the velocity component in z-direction is determined based on the continuity equation, as previously mentioned.
The main advantage of the proposed model is its versatility with respect to the vertical discretization. In this model, vertical spacing does not need to be uniform and can vary according to the problem bathymetry, since the variables u ? (x, y, z, t), v ? (x, y, z, t), and w ? (x, y, z, t) are calculated at any point of the z-axis. The sequence used to calculate three-dimensional flow is summarized as follows:   Table 2 Parameters of the experiment developed by Yu [31] and presented by Jin [11] Case u w (cm/s) u a (m/s) z o (mm)

Applications
An important step in the development of a numerical model consists in the verification of the model against available data sets, which can be obtained from laboratory experiments and from another numerical model already tested. For this model, two test cases were undertaken.

Wind-induced circulation
The three-dimensional model was validated by comparing results obtained by MOTRID with laboratory data sets. An experimental study of turbulent flows induced by wind was conducted by Yu [31] and described by Jin [11]. The test objective is the three-dimensional MOTRID code validation, demonstrating its applicability in wind-forced channel analyses. The numerical velocities profiles obtained by MOTRID are compared with the velocities profiles measured from the experiment by Yu [31], and with profiles obtained by Rosso and Rosman [23], for a three-dimensional model that uses sigma coordinate in the vertical discretization, FIST3D. The experiment by Yu [31] takes into account 24 combinations of air flow and water velocities for different wind conditions. The experiment consists of a channel approximately 0.38 m long with 0.20 m water deep, with a cross section of 0.80 m 9 0.59 m. The wind on the surface generates uniform and constant stress along the entire channel. The air current velocity, measured at 0.10 m above the water surface, varies from 5.7 to 8.0 m/s. A pumping system ensures steady flow in the channel, with cross-sectional averaged velocities varying from 4.75 to 21.4 cm/s. During the experiment, wind direction agrees or not with the flow direction. Figure 5 shows a schematic diagram of the experiment.
All tests were carried out in the hydraulically smooth regime. The measuring section, where profiles were determined, was located near the middle of the channel, where the water depth exhibited a little effect caused by the wind set-up on the free surface. Jin [11] describes only the most significant experiments, listed in Table 2, showing the mean velocity ( u w ), the wind velocity measured 0.10 m above the free surface (u a ), and a roughness parameter z o , depicting equivalent bottom roughness defined as e ¼ 15z o : Numerical simulations were carried out using a grid containing 20 elements, 105 nodes, with constant spacing The Root-Mean-Square-Error (RMSE) values were determined between prediction and measured values. The RMSE error was found to be 10.26%, for case 1, and 8.46%, for case 2. With higher velocities, RMSE decreased: 2.67%, 3.05%, 3.06%, 5.54%, 2.30%, and 2.88% for cases 3, 4, 5, 6, 7, and 8, respectively. In cases 7 and 8, the higher velocities could generate numerical instability in the model. However, the results show that the MOTRID model is accurate in simulating flows with a flat bottom and subject to wind action. The type of flow represented by this experiment is common in regions subject to wind action, easily found in nature.

Channel with non-flat topography
To evaluate the model performance and to demonstrate its three-dimensional predictive capacity, MOTRID was applied to a non-flat topography channel, with two-dimensional behaviour in profile (x-z plane) and no gradient in the y-direction. The channel has approximately 960 m long with 30 m wide.
Numerical simulations were carried out using a grid containing 72 elements, 343 nodes, with constant spacing of 3.53 m in the x-and y-directions. Forty levels were considered in the vertical direction. The channel bathymetry used in the experiment is illustrated in Fig. 14  The results are shown for five different sections, at time equal to 2000s; the sections are shown in Fig. 15. The velocity gradient is maximum in Sect. 2, whereas Sect. 3 presents a slightly lower velocity gradient. The other sections do not present any velocity gradients in the x-direction. Figures 16,17,18,19,20 show a good agreement between MOTRID and FIST3D results; the FIST3D results are assumed as the reference ones. In Sect. 1 (Fig. 16), where there is no variation in bathymetry, the results obtained are very close and depend only on the local variables of each node of the two-dimensional mesh. In Sects. 2 and 3 (Figs. 17 and 18, respectively), the bathymetry variations produce differences in the two formulations, mainly near the bottom. The differences can be attributed to the treatment of the advective terms near the bottom which, in these sections, is significantly different for the two models. Sections 4,5 (Figs. 19 and 20,respectively) show the same differences observed in the previous sections.

Conclusion
This work presents a three-dimensional hydrodynamic model for wind-driven circulation, with a two-dimensional module, which calculates the free surface position and the vertically integrated velocity fields, and a three-dimensional module, which calculates the velocity profiles, using the moving element method.
Coupling of the two modules is performed via the free surface gradient and bottom shear stress. In addition, velocity profiles calculated in the 3D module are adjusted, so that mean vertical velocities are the same in both modules.
In the proposed model, vertical spacing does not need to be uniform and can vary according to the water bodies' bathymetry, diversely from the models that use the sigma coordinates in the vertical discretization. Therefore, the simulation of water bodies with non-flat topography, or with significant variations of the free surface position, can be done in a more efficient way: as there is no need of a fixed division of the water column, deeper regions can have more calculation points than the shallower ones, with a better adjustment in the vertical discretization. This is the main advantage of the model. Future studies can be developed to ensure applications in natural water bodies with more complex geometries.
The moving element method, used for spatial discretization of the three-dimensional module, proved to be efficient, accurate, and stable in simulating flow through the flat-bottomed channel subject to wind action and, also, on a channel with non-flat topography. It is important to mention that, for this type of flow, the correct discretization of advective terms is essential for describing velocity profiles.