Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method

: Though the discontinuous Galerkin method is attracting more and more applications in many ﬁelds due to its local conservation, high-order accuracy, and ﬂexibility for resolving complex geometries, only a few three-dimensional hydrodynamic models based on the discontinuous Galerkin method are present. In this study, a three-dimensional hydrodynamic model with a σ -coordinate system in the vertical direction is developed. This model is discretized in space using the discontinuous Galerkin method and advanced in time with the implicit-explicit Runge–Kutta method. Numerical tests indicate that the developed model is convergent and can obtain better results with a smaller computational time when a higher approximation order is adopted. Other tests with exact solutions also indicate that the developed model can well simulate the vertical circulation under the effect of surface wind stress and the ﬂow under the combined effect of wind stress and Coriolis acceleration terms. The simulation results of tidal ﬂow in part of Bohai Bay, China, indicate that the model can be used for the simulation of tidal wave


Introduction
With the assumption of hydrostatic pressure distribution, the three-dimensional shallow water equations (3D SWEs) can be derived from the Reynolds-averaged threedimensional Navier-Stokes equations. They can predict the vertical distribution of primitive variables and can be used together with advection-diffusion-type equations to simulate the transport and fate of substances, such as temperature, salt, and sediment.
In the last few decades, 3D SWEs have been intensively studied [1], and many hydrodynamic models have been developed and widely used. A brief categorization of these models, based on the numerical methods and the type of grids used in the horizontal direction, is listed in Table 1. For the models at an early stage, the conservative finite difference method on structured grids was usually used. They can be implemented straightforwardly and are computationally efficient, while they approximate the coastlines as staircases and prevent the flexible implementation of variable resolution. Later, unstructured, grid-based hydrodynamic models, such as the finite element and finite volume method, were brought into the main stream, since the coastlines could be accurately represented to give a reasonable simulation at the regional scale [2], and the area of most concern could be simulated using local refinement grids.
In addition to the numerical methods mentioned above, the discontinuous Galerkin (DG) method has received much attention in recent years. It can be viewed as a hybrid between the finite element method and the finite volume method and enjoys most of the strengths of both methods, such as local conservation, compactness, a high order of accuracy, and good application to unstructured meshes. It has been successfully applied in the numerical discretization of two-dimensional shallow water equations [3], two-dimensional non-hydrostatic shallow water equations [4], and Boussinesq-type equations [5].
In the realm of 3D SWEs with DG solutions, Dawson and Aizinger [25] first developed a three-dimensional barotropic model with a z-coordinate system adopted in the vertical direction and presented the corresponding stability analysis [26]. The vertical eddy viscosity terms in the model are explicitly advanced in time. Due to the stiffness of these terms, the time step must be small enough, which limits the applicability of the model. Later, Aizinger et al. [27] made improvements by introducing the baroclinic terms and treating both the vertical eddy viscosity and the vertical diffusion terms implicitly. This model is called UTBEST3D and has been validated with realistic applications. In the same period, Blaise [28] and Comblen [29] developed the three-dimensional baroclinic model named SLIM and adopted it for the simulation of various baroclinic phenomena. Due to the discontinuous nature of the DG method, the water depth between the adjacent elements is generally not continuous such that the lateral boundaries are mismatched, and the post-process measures are needed to smooth the water depth such that the calculation of numerical flux will not be affected. To circumvent the problem of boundary mismatch between adjacent elements, Gandham [30] mapped the time-varying z-coordinate system to a fixed domain by invoking the σ-coordinate transformation and developed a threedimensional barotropic model. However, only one layer is adopted in the vertical direction, which may lead to serious numerical errors when the water depth is large. Conroy and Kubatko [31] developed a three-dimensional barotropic model with a σ-coordinate system adopted in the vertical direction and developed a fast and accurate method for the calculation of depth-averaged velocity. An arbitrary number of layers can be used in the vertical direction. However, non-linear advective terms are not included, and the vertical eddy viscosity terms are explicitly advanced in time. To the authors' knowledge, there is currently no DG method-based three-dimensional hydrodynamic model that includes non-linear advective terms, treats the vertical eddy viscosity terms implicitly, and can use an arbitrary number of vertical layers in the σ-coordinate system simultaneously.
In this paper, a three-dimensional hydrodynamic model that includes nonlinear advective terms and can use an arbitrary number of vertical layers in the σ-coordinate system is developed based on the DG method. The paper is organized as follows. In Section 2, we provide the governing equations of the 3D SWEs in the σ-coordinate system. Then, we present the space discretization and time stepping of the model in Section 3. Numerical experiments are given in Section 4. Finally, conclusions are drawn in Section 5.

Governing Equations
Using the notations given in Table 2, the governing equations of the three-dimensional hydrodynamic model in the physical domain are as follows [13]: ∂u ∂t ∂v ∂t Table 2. Notations for the governing equations of the three-dimensional hydrodynamic model.

Variables Explanation
u(x, y, z, t) velocity of water in x direction v(x, y, z, t) velocity of water in y direction w(x, y, z, t) velocity of water in z direction g acceleration due to gravity θ angle of geographical latitude magnitude of the angular velocity of the Earth f = 2 sin θ Coriolis parameter K m vertical eddy viscosity coefficient K h horizontal eddy viscosity coefficient ρ 0 density τ sx wind stress in x direction τ sy wind stress in y direction η(x, y, t) the surface elevation u η (x, y, t) surface velocity of water in x direction v η (x, y, t) surface velocity of water in y direction b(x, y) bottom elevation bottom velocity of water in y direction z 0 (x, y, t) half height of the bottommost element u c (x, y, t) velocity of water at z 0 in x direction v c (x, y, t) velocity of water at z 0 in y direction The surface boundary conditions for u, v, and w are Likewise, the bottom boundary conditions for u, v, and w are where C f is determined by matching a logarithmic bottom layer to the model at height z 0 and calculated as [32] To accurately represent the bottom and surface geometry, the σ transformation developed by Philips [33] is adapted According to the principle of chain differentiation, partial derivatives of the vari- ∂U ∂t with where ω is related to w by The surface boundary condition in the computational domain is given by and the bottom boundary condition is given by Integrating Equation (11) from bottom to surface and applying the boundary condition about ω leads to the primitive continuity equation as ∂D ∂t

Numerical Implementation
In this section, we will detail the numerical implementation of our model based on the DG method. The domain partition and the discontinuous polynomial space are first described. Then, the numerical discretization of each part is presented, followed by the time-stepping scheme.

Domain Partition and Polynomial Space
For an arbitrary three-dimensional domain Ω 3d , its horizontal projection is denoted as Ω 2d and is partitioned by Ele2d non-overlapping triangular or quadrilateral elements; i.e., Ω 2d,h = {Ω e }, e ∈ [1, Ele2d]. In the vertical direction, Ω 2d,h is extruded by N L layers to obtain the three-dimensional computational domain Ω 3d,h , and the computational domain consists of Ele2d × N L triangular prisms or quadrangular prisms. In this study, the threedimensional element obtained by the extrusion of Ω e is denoted as Ω (e,L) , L ∈ [1, N L ], with the bottommost element denoted as Ω (e,1) and the upmost element denoted as Ω (e,N L ) . Figure 1 gives a schematic view of two columns of three-dimensional computational elements with Figure 1a showing the quadrangular prisms and Figure 1b showing the triangular prisms. To illustrate the discontinuous nature of the DG method, the discontinuities between the adjacent computational elements have been exaggerated, and no gap exists in reality. In addition, we introduce the following finite dimensional space of polynomials:   In addition, we introduce the following finite dimensional space of polynomials: where L 2 (Ω) is the space of the square-integrable functions and P (N h ,N v ) Ω (e,L) is the complete space of polynomials defined over Ω (e,L) , which is of order at most N h in the horizontal direction and at most N v in the vertical direction. In this study, V h,3d is spanned by N p Lagrangian functions, l i (x), i ∈ 1, N p , and the dimension, N p , depends on the approximation order and element type. For triangular prisms, while for quadrilateral prisms, For any two adjacent elements Ω e1 and Ω e2 , let ε = ∂Ω e1 ∩ ∂Ω e2 be the unique interior interface between these two elements, and let n − and n + be the outward unit normal for these two elements on ε. For scalar field c ∈ V h,3d , let c ± ε be the trace of c on interface ε from the interior of Ω e1 (2) . We further define the mean value {c} and the jump value [c] on interface ε as On the boundary of the computational domain, they are defined as where c is the single value and n is the outward unit normal on the boundary edges. We note that the jump of scalar c is a vector, and we further denote its horizontal and vertical components by [c] h and [c] σ , respectively.

Numerical Discretization of Momentum Equations
For Equation (12), both the Coriolis acceleration terms and the vertical eddy viscosity terms are treated implicitly, while the others are treated explicitly.

Convective and Bottom Topography Terms
For the numerical discretization of Equation (12), we consider the following general formulation: find the local approximate solution vector U h ∈ (V h,3d ) 2 such that for all test functions φ ∈ V h,3d and Ω (e,L) ∈ Ω 3d,h we have where U − h and U + h are the values of the local approximate solution vector U h on the interface of two adjacent elements, and and n · F * U − h , U + h are the flux term and the numerical flux at the element boundary, respectively. For the three-dimensional computational mesh used in this study, the third component of the unit outward normal is zero, i.e., n σ = 0, at the lateral boundary of the computational elements, the calculation of the numerical flux degenerates to that of the two-dimensional shallow water equations, and the Harten-Lax-van Leer with contact (HLLC) Riemann solver [34] is adopted. As n x = n y = 0 at the top and bottom boundaries of the computational elements, only the vertical convection needs to be considered, and the upwind flux is taken, i.e., where sgn is the sign function. At the top boundary of the computational domain (σ = 0) and bottom boundary of the computational domain (σ = −1), the numerical flux is set to zero. (25) are the discretization corresponding to the horizontal eddy viscosity terms and the vertical and Coriolis acceleration terms, which will be detailed in the forthcoming parts. The requirement that Equation (25) is valid for all φ ∈ V h,3d and can be imposed by considering Equation (25) for all the basis functions that span the functional space V h,3d , i.e., setting For the discretization of both the volume and the surface terms, the quadrature-free manner of Hesthaven and Warburton [35] is adopted. We assume that these terms also lie in the finite dimensional space (V h,3d ) 2 and can be expressed as where l f l i l f dx,

Horizontal Eddy Viscosity Terms
For the discretization of the horizontal eddy viscosity terms in Equation (12), the local discontinuous Galerkin method by Cockburn and Shu [36] is adopted. To discretize the horizontal eddy viscosity terms in the momentum equation about Du, we first express it in the following form: where r is an auxiliary function. Following the philosophy of the local discontinuous Galerkin method, we introduce the auxiliary variables and cast Equation (30) into the following form:q = q x ,q y = ∇ xy u, with υ h = K h D and ∇ xy = (∂/∂x, ∂/∂y). For the numerical discretization of Equation (31), where the terms marked with superscript "*" are the flux terms and are given as where τ is the penalty parameter defined as where N = max(N h , N v ), n 0 is the number of neighbors of the element, i.e., five for triangular prisms and six for quadrilateral prisms, A is the area of the interface, and V is the average volume of the two adjacent elements.
and express these terms aŝ where n h = n x , n y is the vector of the horizontal components of the outward unit normal outward at the facial interpolation points. The right-hand side of Equation (38) corresponds to the numerical discretization of the second-order term for the horizontal momentum equation about Du, and the same procedure applies to that about Dv.

Vertical Eddy Viscosity and Coriolis Acceleration Terms
For the discretization of the vertical eddy viscosity and Coriolis acceleration terms in Equation (12), we first express it in the following form: where g Du N is the Neumann boundary condition about Du on Ω N and is given in Equations (17) and (18); g Du D is the Dirichlet boundary condition about Du on Ω D . The symmetric interior penalty discontinuous Galerkin (SIPG) method [37] is adopted, and the corresponding primal form is given as where υ v = K m /D 2 is the vertical eddy viscosity coefficient in the computational domain.
The primal form about Dv is similar to that about Du. In this study, the Neumann boundary conditions for Du and Dv given in Equation (18) are treated implicitly, and the bottom velocity module u 2 c + v 2 c is linearized by using the old values to simplify the calculation. As Equation (40) is only solved in the vertical direction, this equation is independent for each column of prisms and can be solved independently.

Numerical Discretization of the Primitive Continuity Equation
The water depth is calculated according to Equation (19), and the semi-discrete form is given as where M e , S e = S x,e , S y,e and M ∂ f e are the mass matrix, stiff matrix and edge-mass matrix defined over the two-dimensional computational element Ω e , respectively. Their definitions are analogous to the three-dimensional counterpart defined in Equation (29). two-dimensional computational element Ω e and are obtained through the depth integration of their three-dimensional counterpart to preserve the local mass conservation properties.

Calculation of Vertical Velocity
According to Equations (11) and (19), the governing equation for vertical velocity ω is given as ∂ω ∂σ For the space discretization of Equation (42), suppose that for all φ ∈ V h,3d and Ω (e,L) ∈ Ω h,3d , we have where a h is the approximate solution to a and belongs to V h,3d . In Equation (43), the boundary of Ω (e,L) is split into three parts, i.e., the lateral boundary ∂Ω Lat (e,L) , the top boundary ∂Ω Top (e,L) , and the bottom boundary ∂Ω Bot (e,L) , and terms marked with an asterisk are the numerical flux terms. Here, the numerical flux n is calculated according to the HLLC Riemann solver, and n x (DU) * h + n y (DV) * h is the depth-averaged counterpart. At the bottom boundary, ∂Ω Bot (e,L) , ω * h is taken as follows: The numerical flux is taken to be the boundary condition at the bottom in the case of the bottommost element or using the solution from the element below. Thus, the integral defined over ∂Ω Top (e,L) diminishes. With such a definition of numerical flux about ω * h given in Equation (44), we can obtain the vertical velocity ω h layer by layer in each water column starting at the bottom and using the solution from the element below to provide an initial condition.

Time Stepping
In the discretization presented above, the ordinary differential equations Equation (28) and Equation (41) are obtained for the momentum equations and the primitive continuity equation, and they can be formulated as where y(t) denotes the vector of all discrete degrees of freedom of a step and f exp l (y(t), t) and f impl (y(t), t) denote the terms treated explicitly and implicitly. For the momentum equations, the explicit terms consist of nonlinear advection, horizontal eddy viscosity, and bottom topography source terms, while the implicit terms consist of the vertical eddy viscosity and Coriolis acceleration terms. For the primitive continuity equation, all terms are treated explicitly. In this study, the implicit-explicit Runge-Kutta method is adapted to integrate system (45) in time. For an s-stage implicit Runge-Kutta method defined by coefficients a i,j , b i , and c i and a σ = s + 1 stage explicit Runge-Kutta method defined by coefficientsâ i,j ,b i , andĉ i , this time-stepping method is formulated as where y n are the known values at time t n , ∆t is the time step, , and y n+1 are the unknown values at t n+1 = t n + ∆t. The time-stepping method adopted is derived by Ascher [38], and the corresponding coefficients are The time-stepping method adopted is derived by Ascher [38], and the corresponding coefficients are 1.
Calculate the explicit terms

Numerical Tests
In this section, several numerical experiments are conducted to verify the performance of the developed model. For all the experiments, the acceleration due to

Numerical Tests
In this section, several numerical experiments are conducted to verify the with γ = 2 − √ 2 /2 and δ = 1 − 1/(2γ). The global time-stepping algorithm from t n to t n+1 = t n + ∆t is as follows: 1.
Calculate the explicit terms K Calculate the intermediate water depth as Likewise, the intermediate conservative variables Calculate the water depth D n+1 and the conservative variables U n+1 at t n+1 as

4.
Integrate U n+1 along the water depth to obtain the final depth-averaged momenta (DU) n+1 and (DV) n+1 , followed by the calculation of the final vertical velocity ω n+1 .

Numerical Tests
In this section, several numerical experiments are conducted to verify the performance of the developed model. For all the experiments, the acceleration due to gravity, g, is set to 9.81 m/s 2 . In addition, all the experiments are run on an Intel Xeon E5-2620 processor with 16 GB of internal memory. Our program is parallelized using OpenMP to run on six cores.

Manufactured Solution
This test is used to verify the convergence property of the developed model. As the 3D SWEs are a non-linear system, designing a nontrivial function satisfying the equations and the boundary conditions at the same time is not a trivial task. Following the philosophy presented by Salari and Knupp [39], we manually give the solution that satisfies the continuity equation and solve a Dirichlet problem for a modified system that includes a forcing term.

(52)
The vertical velocity at an arbitrary location is obtained through the integration of the continuity equation, Equation (11), from the bottom to the studied location. The horizontal projection of the horizontal domain is Ω 2d = [0, 100] × [0, 100], and the whole simulation lasts for 86.4 s.
Tests are run with successively refined meshes. For the coarsest mesh, there is one element in the horizontal direction and one layer in the vertical direction. We refine the mesh up to four times by partitioning each quadrilateral prism into eight and compute the L 2 error of the numerical results against the analytical solutions as where Q S e (x) is the numerical result and Q exact e (x) is the analytical solution. To measure the convergence properties of the model, we further define the convergence rate (CR) as follows where L m 2 and L m−1 2 are the numerical errors on two successive refined meshes and L m and L m−1 are the characteristic lengths of the meshes, which are set to be the longest edge length of the quadrilateral prisms. Tables 3 and 4 give the L 2 error for each variable and the corresponding CR for N h = N v = 1 and 2, respectively. N e denotes the number of elements in the horizontal direction and N L the number of vertical layers. We note that for water depth D, horizontal momenta Du and Dv converge at the optimal rate, i.e., CR = 2 for the former case and 3 for the latter, whereas the vertical velocity ω converges at the suboptimal rate for both cases. We owe this to the fact that the vertical velocity is computed from the numerical results for Du, Dv, and their depth-averaged counterparts using the continuity Equation (42), and the CR for vertical velocity would be influenced as the horizontal momenta are not exact, as pointed out by Dawson and Aizinger [25]. Table 3. N h = N v = 1, the L 2 error for each variable, and the corresponding CR.  Figure 2 shows the relation between the L 2 error and DOFs and the normalized CPU time. Figure 2a indicates that the developed model is convergent, as the L 2 error for each variable diminishes in both cases as the computational mesh is refined. Figure 2b indicates that we can achieve a given error tolerance with fewer DOFs on a coarser mesh when a higher approximation order is adapted. For example, when the error tolerance about water depth D is taken as 1.0 × 10 −6 , the normalized CPU time is about 0.033 when the approximation order is N h = N v = 2, whereas it is about 0.18 when the approximation order is N h = N v = 1.

Tide Wave Propagation in a Semi-Closed Bay
This test is used to indicate the accuracy of the developed model for the sim of both the horizontal and vertical velocity fields. Figure 3 shows the horizontal l of Gauge points A (2500, 5000), B (52,500, 5000), and C (92,500, 5000) and a sketc horizontal domain. This domain is 100,000 m long and 10,000 m wide and is par

Tide Wave Propagation in a Semi-Closed Bay
This test is used to indicate the accuracy of the developed model for the simulation of both the horizontal and vertical velocity fields. Figure 3 shows the horizontal locations of Gauge points A (2500, 5000), B (52,500, 5000), and C (92,500, 5000) and a sketch of the horizontal domain. This domain is 100,000 m long and 10,000 m wide and is partitioned by a uniform quadrilateral element with a size of 1000 m in both directions. The still water depth is h = 12 m. A cosine tidal wave with amplitude A = 0.25 m and period T = 12 h is imposed at the left boundary, i.e., x = 0, whereas the other boundaries are all treated as slip walls. Both the eddy viscosity terms and the Coriolis acceleration terms are not included.
is the frequency. For the initial condition, all the velocity compone are set to zero, and the surface elevation is set according to the analytical solution. T number of vertical layers is set to 12, and we run this simulation for 24 h. Figure 4 gives the comparison between the numerical result and the analyti solution for surface elevation , horizontal velocity , and vertical velocity at Gaug A, B, and C, and the vertical coordinate for the three Gauges is =−2.5 m. It can be se that the numerical solution is generally in good agreement with the analytical soluti indicating that the three-dimensional model developed in this study has good accuracy the simulation of water level and velocity. According to Neuman and Pierson [40], the analytical solution of the surface elevation η, the horizontal velocities u and v, and the vertical velocity w at a point P(x, y, z) at time t are given as where ω = 2π/T is the frequency. For the initial condition, all the velocity components are set to zero, and the surface elevation is set according to the analytical solution. The number of vertical layers is set to 12, and we run this simulation for 24 h. Figure 4 gives the comparison between the numerical result and the analytical solution for surface elevation η, horizontal velocity u, and vertical velocity w at Gauges A, B, and C, and the vertical coordinate for the three Gauges is z = −2.5 m. It can be seen that the numerical solution is generally in good agreement with the analytical solution, indicating that the three-dimensional model developed in this study has good accuracy in the simulation of water level and velocity. Figure 4 gives the comparison between the numerical result and the analytical solution for surface elevation , horizontal velocity , and vertical velocity at Gauges A, B, and C, and the vertical coordinate for the three Gauges is =−2.5 m. It can be seen that the numerical solution is generally in good agreement with the analytical solution, indicating that the three-dimensional model developed in this study has good accuracy in the simulation of water level and velocity.

Wind-Induced Water Circulation
In this test, the wind-driven circulation in a rectangular closed basin is simulated. Following Liu et al. [41], we run this case to test the accuracy of the proposed numerical model in predicting the vertical stratified circulation.
With a homogeneous Dirichlet boundary condition adapted at the bottom boundary, Koutitas et al. [42] derived the analytical solution of the horizontal velocity in the center of the basin as where ℎ is the still water depth of the basin. Likewise, Huang [43] also derived that when the Neumann boundary condition is adapted, the result reads where | | is the velocity magnitude at the center of the bottommost element, and | | = 0.005 / in this study. In this part, these two situations are simulated.
This case is set as follows: the basin is 2000 m long and 800 m wide, ℎ = 10 m, and

Wind-Induced Water Circulation
In this test, the wind-driven circulation in a rectangular closed basin is simulated. Following Liu et al. [41], we run this case to test the accuracy of the proposed numerical model in predicting the vertical stratified circulation.
With a homogeneous Dirichlet boundary condition adapted at the bottom boundary, Koutitas et al. [42] derived the analytical solution of the horizontal velocity in the center of the basin as where h is the still water depth of the basin. Likewise, Huang [43] also derived that when the Neumann boundary condition is adapted, the result reads where |u c | is the velocity magnitude at the center of the bottommost element, and C f |u c | = 0.005 m/s in this study. In this part, these two situations are simulated.
This case is set as follows: the basin is 2000 m long and 800 m wide, h = 10 m, and τ sx = 1.5 Pa. The Coriolis acceleration and horizontal eddy viscosity terms are ignored, and the vertical eddy viscosity coefficient K m is set to be 0.001 m 2 s −1 ; the whole simulation lasts for 3600 s. Figure 5a,b shows the comparison between the numerical results and the analytical solution for the vertical distribution of horizontal velocity at the center of the basin. The comparison indicates that the numerical results match well with the analytical solution for both cases.  Table 5 gives the RMSE for diffe simulations. We find that the more vertical layers there are, the better the result will b each case.  We also run the simulations with different numbers of vertical layers (N L ). To quantitatively measure the difference between different simulations, we define the root mean square error (RMSE) as where N T is the number of simulated results, u j is the numerical result at point j, and u a,j is the analytical solution at the same point. Table 5 gives the RMSE for different simulations. We find that the more vertical layers there are, the better the result will be for each case.

Generation of the Ekman Profile
Following Trahan et al. [44], we run this test to check whether the developed model can generate the Ekman layer, where the Coriolis acceleration terms, the pressure gradient, and the turbulent drag are in balance. The analytic solution for the Ekman velocity profile was developed by Ekman [45] and is expressed as where is the folding depth, = 0.1 m 2 /s, = 0.1 Pa, and is the latitude ( 45°). Theoretically, the flow direction at the surface w counterclockwise 45 degrees to that of the wind stress.
The computational domain is 40,000 km long and 40,000 km wide and by a uniform quadrilateral element with a size of 2000 km in both direction depth of 200 m is set, and 50 vertical layers are used. The bottom drag coeffic to 0.0025. We run this simulation for 200 days. Figure 6 gives the comp analytic solution and the computed Ekman profile. The numerical resu match the analytic solution and individual velocity components with an av less than 0.001 m/s.

Tidal Flow in Bohai Bay
For the last case, we simulate tidal flow in the western region of Bohai B ability of the developed model for the simulation of actual tidal wave moti The Bohai Bay is one of the three major bays of the Bohai Sea in China. It is west of the Bohai Sea and connected with Hebei and Shandong Province and It's a typical semi-closed sea area. The average tidal range in Bohai Bay is tidal current is irregularly semi-diurnal. In recent years, a lot of port and projects in Bohai Bay have been carried out, which has led to a large number engineering structures. The accurate simulation of the tidal current fiel significance for the convection-diffusion process of temperature, salt, pollu well as the sediment transport in Bohai Bay.
The domain geometry for this case is shown in Figure 7a; the still water

Tidal Flow in Bohai Bay
For the last case, we simulate tidal flow in the western region of Bohai Bay to test the ability of the developed model for the simulation of actual tidal wave motion problems. The Bohai Bay is one of the three major bays of the Bohai Sea in China. It is located in the west of the Bohai Sea and connected with Hebei and Shandong Province and Tianjin City. It's a typical semi-closed sea area. The average tidal range in Bohai Bay is 2.5 m, and its tidal current is irregularly semi-diurnal. In recent years, a lot of port and reclamation projects in Bohai Bay have been carried out, which has led to a large number of large-scale engineering structures. The accurate simulation of the tidal current field is of great significance for the convection-diffusion process of temperature, salt, pollutants, etc., as well as the sediment transport in Bohai Bay.
The domain geometry for this case is shown in Figure 7a; the still water depth ranges from 38 m at the east open boundary to 2 m at the coast. There are five measuring gauges in the computational domain, and the distribution of these points is shown in Figure 7b. The computational domain is partitioned by 14,627 triangle grids, and the edge length of the finite element mesh ranges from 2 km to 100 m. Five equidistant layers are used in the vertical direction. Wetting/drying is not accounted in the present study.

of 24
in the computational domain, and the distribution of these points is shown in Figure 7b. The computational domain is partitioned by 14,627 triangle grids, and the edge length of the finite element mesh ranges from 2 km to 100 m. Five equidistant layers are used in the vertical direction. Wetting/drying is not accounted in the present study. The tidal elevations at the east open boundaries were provided by Chinatide [46], which is a tidal forecast system through nine harmonic constants of Q1, P1, O1, K1, N2, M2, S2, K2, and Sa. The Coriolis parameter was taken as 9.1557 × 10 s , and the horizontal eddy viscosity coefficient was set to 70 m 2 s −1 . The vertical eddy viscosity coefficient was modelled using the − ϵ turbulence closures provided by GOTM [47], and the coupling between GOTM and the hydrodynamic model is established following the same philosophy presented by Tuomas et al. [32]. The bottom roughness length is set as 0.1 mm. Field measurements from 12/08/2018 to 12/09/2018 are used to verify the developed model. Figure 8 shows the flow field at the surface and the middle and bottom layers at the time of flood tide (left) and ebb tide (right), respectively. As we can see, there is little difference between the velocity fields of the surface and the middle layers. Because of the bottom friction, the velocity of the bottom layer is smaller than that of the surface and the middle layers, which is about 60% of the magnitude of the latter two. The results of the three-dimensional velocity field show that the model can reflect the three-dimensional characteristics of tidal flow. The tidal elevations at the east open boundaries were provided by Chinatide [46], which is a tidal forecast system through nine harmonic constants of Q 1 , P 1 , O 1 , K 1 , N 2 , M 2 , S 2 , K 2 , and Sa. The Coriolis parameter was taken as 9.1557 × 10 −5 s −1 , and the horizontal eddy viscosity coefficient was set to 70 m 2 s −1 . The vertical eddy viscosity coefficient was modelled using the k − turbulence closures provided by GOTM [47], and the coupling between GOTM and the hydrodynamic model is established following the same philosophy presented by Tuomas et al. [32]. The bottom roughness length is set as 0.1 mm. Field measurements from 12/08/2018 to 12/09/2018 are used to verify the developed model. Figure 8 shows the flow field at the surface and the middle and bottom layers at the time of flood tide (left) and ebb tide (right), respectively. As we can see, there is little difference between the velocity fields of the surface and the middle layers. Because of the bottom friction, the velocity of the bottom layer is smaller than that of the surface and the middle layers, which is about 60% of the magnitude of the latter two. The results of the three-dimensional velocity field show that the model can reflect the three-dimensional characteristics of tidal flow. Figure 9 shows the comparison between the simulated surfaced elevation and the field measurements at Gauges P1 and P2. The simulated data are quantitatively very close to the field measurements, and the tidal phasing is virtually identical. The left column of Figure 10 gives the comparison between the simulated velocity and the field measurements at gauges P3, P4, and P5. Small discrepancies exist between the simulated velocity and the field measurements. This is especially the case for P5, where the maximum simulated velocity is about 15% smaller than that of the field measurements. It can be attributed to that the computational mesh used in this simulation is relatively coarse, and the wetting/drying is not considered. Although differences exist, the simulated data are in phase with the field measurements as the comparison between the simulated flow direction and the field measurements shows in the right column of Figure 10. Generally, the developed hydrodynamic model can be used for the simulation of tidal wave motion in realistic situations.  Figure 9 shows the comparison between the simulated surfaced elevation and the field measurements at Gauges P1 and P2. The simulated data are quantitatively very close to the field measurements, and the tidal phasing is virtually identical. The left column of Figure 10 gives the comparison between the simulated velocity and the field maximum simulated velocity is about 15% smaller than that of the field measurements. It can be attributed to that the computational mesh used in this simulation is relatively coarse, and the wetting/drying is not considered. Although differences exist, the simulated data are in phase with the field measurements as the comparison between the simulated flow direction and the field measurements shows in the right column of Figure 10. Generally, the developed hydrodynamic model can be used for the simulation of tidal wave motion in realistic situations. maximum simulated velocity is about 15% smaller than that of the field measurements. It can be attributed to that the computational mesh used in this simulation is relatively coarse, and the wetting/drying is not considered. Although differences exist, the simulated data are in phase with the field measurements as the comparison between the simulated flow direction and the field measurements shows in the right column of Figure 10. Generally, the developed hydrodynamic model can be used for the simulation of tidal wave motion in realistic situations.

Conclusions
In this study, a three-dimensional hydrodynamic model based on the Discontinuous Galerkin method is developed. The σ-coordinate system is used in the vertical direction to circumvent the mismatch of lateral faces due to the discontinuous nature of the solution between different elements. Non-linear advective terms are included, and the model is marched in time with the implicit-explicit Runge-Kutta method. The vertical eddy viscosity and the Coriolis acceleration terms are discretized with the symmetric interior penalty discontinuous Galerkin method and are treated implicitly, while others are treated explicitly.
The developed model is validated with a series of numerical tests. The case with the manufactured solution indicates that the model is convergent, and the L 2 error about the variables diminishes as we increase the approximation order or refine the computational mesh. Calculating the CR indicates that water depth D and horizontal momenta Du and Dv can all converge at the optimal order, whereas the vertical velocity converges in a suboptimal order. This case also indicates that we can obtain a better numerical result on a computational mesh with fewer DOFs when a higher approximation is adapted. Cases of tide-induced three-dimensional flow in a semi-closed bay, wind-induced water circulation, and the generation of the Ekman profile indicate that the developed model can well simulate the velocity field in both horizontal and vertical directions, can reasonably simulate the vertical stratified circulation, and can simulate water motion under the combined effect of wind stress and Coriolis acceleration terms. For the Bohai Bay case, the simulated surface elevation, velocity, and flow direction match reasonably well with the field measurements, indicating the model can be used for the simulation of tidal wave motion in realistic situations.
More applications of the developed model and development of the wetting/drying treatment and the baroclinic model will be presented in forthcoming work.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.