A new three-dimensional terrain-following tidal model of free-surface flows

Abstract

A three-dimensional hydrodynamic model is presented which combines a terrain-following vertical coordinate with a horizontally orthogonal curvilinear coordinate system to fit the complex bottom topography and coastlines near estuaries, continental shelves, and harbors. To solve the governing equations more efficiently, we improve the alternating direction implicit method, which is extensively used in the numerical modeling of horizontal two-dimensional shallow water equations, and extend it to a three-dimensional tidal model with relatively little computational effort. Through several test cases and realistic applications, as presented in the paper, it can be demonstrated that the model is capable of simulating the periodic to-and-fro currents, wind-driven flow, Ekman spirals, and tidal currents in the near-shore region.

Keywords

three-dimensional numerical model / alternating direction implicit method / tidal flow / orthogonal curvilinear coordinates / terrain-following

Introduction

Numerical models are widely used at present to simulate free surface flows in rivers, estuaries, and oceans. With the rapid increase in computer power in recent years, the use of three-dimensional (3D) hydrodynamic models has shifted from academic research to practical applications. Such a model supplies the basic current structure for pollutant transport, oil-spill movement, and suspended sediment transport models ( Drago et al., 2000; Hu et al., 2009).

Compared with the depth-averaged two-dimensional (2D) model, which has been extensively used for simulating storm surges, pollutant transport, and the evolution of rivers ( Chen et al., 2009; Zhang et al., 2010; Hu et al., 2011b), a 3D model is better able to reflect vertical variations in velocity and allows the bed friction stress to be determined more accurately ( Mao et al., 2008; Xie and Zhuang, 2010; Zhao et al., 2012). Therefore, the tidal field obtained using a 3D model is typically more realistic and reasonable than that obtained using a 2D approach.

In the context of 3D simulation, it has often been noted that the bottom topography of the sea is irregular and that the ordinary z-coordinate system encounters certain disadvantages in following it. A staircase representation of bathymetry leads to numerical errors and may not accurately represent bottom boundary layer flows. Gerdes ( 1993) found that neither barotropic nor baroclinic coastal-trapped waves could be accurately represented using a z-coordinate model. Bell ( 1999) also demonstrated that z-coordinate models can generate significant vorticity errors. Therefore, terrain-following coordinates are widely used in tidal modeling. One terrain-following coordinate, the σ coordinate, was first introduced for use in atmospheric modeling ( Phillips, 1957) and has become a standard vertical coordinate for the modeling of tides or other patterns of coastal flows ( Zhang et al., 2006; Young et al., 2007; Keilegavlen and Berntsen, 2009; Hu et al., 2011a). The main advantage of the σ coordinate is that when it is cast in a different form, a smooth representation of the bottom topography is obtained; one can also easily incorporate both a bottom boundary layer and a surface boundary layer into the coordinate system.

In the horizontal direction, traditional Cartesian coordinate grids are not well suited for application to a coastal region with complex coastlines and bathymetry; such grids can lead to poor accuracy of the numerical solution and a non-optimal compromise between accuracy and computer resources. By comparison, orthogonal curvilinear coordinates can accommodate a complex coastline more precisely. Through a coordinate transformation, such an irregular physical grid can be transformed into a regular computational one, and the grid can fit any boundary shape of flowing fields, thereby improving the calculation precision at the boundary with only a slight sacrifice in computational time as a result of the need for coordinate transformation ( Ly et al., 1999; Hu et al., 2009). Indeed, this transformation effort is minimal compared with that required for an unstructured grid, which is also known for its adaptability.

At present, most existing 3D models assume similar simplifications to the Reynolds-averaged Navier-Stokes (RANS) equations, based on the Boussinesq approximation and hydrostatic pressure assumption, and differ in their numerical schemes. For instance, the governing equations can be discretized using various methods such as the finite difference, finite element, or finite volume method. The time integration schemes range from fully explicit to fully implicit. Currently, most existing 3D models use semi-implicit schemes because they use an explicit scheme for the horizontal direction but an implicit scheme for the vertical direction ( Blumberg and Mellor, 1987; Toro and Gomez, 1998). In recent years, several fully implicit splitting methods have been developed and applied ( Casulli and Stelling, 1998; Stansby and Zhou, 1998).

In terms of time integration, many of the primitive equation models in widespread use are solved explicitly. To overcome the severe time-step restriction imposed by the Courant−Friedrichs−Lewy (CFL) criterion on the solution of the water elevation of the free surface, a mode-splitting technique is used ( Lazure and Dumas, 2008). This technique separates the entire system into external and internal modes and solves each of them separately at different time steps; thus, a small time step can be used on the barotropic equations, and a larger time step can be used on the baroclinic equations. This method has the disadvantages that terms must be introduced to account for the interaction between the two modes and that consistency between the internal and external modes must be maintained. Many models, such as the Princeton Ocean Model (POM), S-coordinate Primitive Equation Model (SPEM) and Miami Isopycnic Coordinate Ocean Model (MICOM), use this technique. By contrast, the tidal model described here solves the continuity and momentum equations using the semi-implicit alternating direction implicit (ADI) method, which was first applied to simulate the horizontal tidal current field by Leendertse (1970) and has been extended to 3D models in the past decade ( Yang and Ozer, 1997; Shi et al., 1999). The advantages of the ADI method are its stability and high efficiency, which mean that typically Courant numbers of 5‒10 times larger can be used for the barotropic modes. In this paper, a fairly simple method is used to extend the ADI method from 2D to 3D. It is demonstrated that the ADI method can be used to solve for 3D coastal flows with good stability.

The paper is organized as follows: In section 2, the model and its boundary conditions are introduced, and in section 3, the numerical scheme is proposed. Subsequently, the performance of the model is tested in section 4, including both ideal testing and practical application. Finally, conclusions are offered in section 5.

The governing equations

Hydrodynamic equations

The equations that form the basis of the tidal model describe the velocity and surface elevation field. Consider the hydrodynamic systems for free surface flows given orthogonal curvilinear coordinates in the horizontal direction and the σ coordinate in the vertical direction, where the bottom is at σ=‒1 and the surface is at σ=0(see Fig. 1). The processes under study occur at a horizontal scale that is at least one order of magnitude larger than the vertical scale; thus, the hydrostatic assumption is used to simplify the equations.

The continuity equation is

(1) h_{\mathrm{1}}{h}_{\mathrm{2}}{\displaystyle \frac{\partial \varsigma }{\partial t}}+{\displaystyle \frac{\partial (h_{\mathrm{2}}uD)}{\partial \epsilon }}+{\displaystyle \frac{\partial (h_{\mathrm{1}}vD)}{\partial \eta }}+h_{\mathrm{1}}{h}_{\mathrm{2}}{\displaystyle \frac{\partial w}{\partial \sigma }}=\mathrm{0.}

The ε-momentum equation is

(2) {\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{u}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial u}{\partial \epsilon }}+{\displaystyle \frac{v}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial u}{\partial \eta }}+{\displaystyle \frac{\omega }{D}}{\displaystyle \frac{\partial u}{\partial \sigma }}-\tilde{f}v-fv=-{\displaystyle \frac{g}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial \varsigma }{\partial \epsilon }}+A_H\left({\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial {\epsilon }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial {\eta }^{\mathrm{2}}}}\right)+{\displaystyle \frac{\mathrm{1}}{D^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \sigma }}\left(Av{\displaystyle \frac{\partial u}{\partial \sigma }}\right)

The η-momentum equation is

(3) {\displaystyle \frac{\partial v}{\partial t}}+{\displaystyle \frac{u}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial v}{\partial \epsilon }}+{\displaystyle \frac{v}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial v}{\partial \eta }}+{\displaystyle \frac{\omega }{D}}{\displaystyle \frac{\partial v}{\partial \sigma }}+\tilde{f}u+fu=-{\displaystyle \frac{g}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial \varsigma }{\partial \eta }}+A_H\left({\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial {\epsilon }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial {\eta }^{\mathrm{2}}}}\right)+{\displaystyle \frac{\mathrm{1}}{D^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \sigma }}\left(Av{\displaystyle \frac{\partial v}{\partial \sigma }}\right)

where u, v, and w are the ε, η, and σ components of the current; \varsigma is the free surface elevation; f is the Coriolis parameter; A_H and Av are the horizontal and vertical eddy viscosity coefficients; h_{\mathrm{1}} and h_{\mathrm{2}} are the Lame coefficients along the ε and η directions, ( Blumberg and Herring, 1987); g is the gravitational acceleration; D=\varsigma +h is the total water depth; and h denotes the bottom bathymetry. \tilde{f}={\displaystyle \frac{v}{h_{\mathrm{1}}{h}_{\mathrm{2}}}}{\displaystyle \frac{\partial h_{\mathrm{2}}}{\partial \epsilon }}-{\displaystyle \frac{u}{h_{\mathrm{1}}{h}_{\mathrm{2}}}}{\displaystyle \frac{\partial h_{\mathrm{1}}}{\partial \eta }} denotes the curvature term.

Integrating Eqs. (1), (2), and (3) over the depth and applying a kinematics condition at the free surface leads to the following equations:

(4) h_{\mathrm{1}}{h}_{\mathrm{2}}{\displaystyle \frac{\partial \varsigma }{\partial t}}+{\displaystyle \frac{\partial (h_{\mathrm{2}}UD)}{\partial \epsilon }}+{\displaystyle \frac{\partial (h_{\mathrm{1}}VD)}{\partial \eta }}=\mathrm{0},

(5) {\displaystyle \frac{\partial U}{\partial t}}=-{\displaystyle \frac{g}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial \varsigma }{\partial \epsilon }}+{\displaystyle \frac{{\tau }_{s\epsilon }}{\rho D}}-{\displaystyle \frac{{\tau }_{b\epsilon }}{\rho D}}+{\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{0}}{M}_{\epsilon }d\sigma },

(6) {\displaystyle \frac{\partial V}{\partial t}}=-{\displaystyle \frac{g}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial \varsigma }{\partial \eta }}+{\displaystyle \frac{{\tau }_{s\eta }}{\rho D}}-{\displaystyle \frac{{\tau }_{b\eta }}{\rho D}}+{\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{0}}{M}_{\eta }d\sigma },

where the vertically integrated horizontal velocities are defined as

(7) (U,V)={\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{0}}(u,v)d\sigma },

and

(8) M_{\epsilon }=-\left({\displaystyle \frac{u}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial u}{\partial \epsilon }}+{\displaystyle \frac{v}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial u}{\partial \eta }}+{\displaystyle \frac{\omega }{D}}{\displaystyle \frac{\partial u}{\partial \sigma }}\right)+A_H\left({\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial {\epsilon }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}u}{\partial {\eta }^{\mathrm{2}}}}\right)+fv+\tilde{f}v,

(9) M_{\eta }=-\left({\displaystyle \frac{u}{h_{\mathrm{1}}}}{\displaystyle \frac{\partial v}{\partial \epsilon }}+{\displaystyle \frac{v}{h_{\mathrm{2}}}}{\displaystyle \frac{\partial v}{\partial \eta }}+{\displaystyle \frac{\omega }{D}}{\displaystyle \frac{\partial v}{\partial \sigma }}\right)+A_H\left({\displaystyle \frac{\mathrm{1}}{h_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial {\epsilon }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{h_{\mathrm{2}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}v}{\partial {\eta }^{\mathrm{2}}}}\right)-fu-\tilde{f}u.

{\tau }_{s\epsilon } and {\tau }_{s\eta } are the wind stress components, and {\tau }_{b\epsilon } and {\tau }_{b\eta } are the bottom stress components, which will be given in section 2.2.

Subtracting Eq. (5) from Eq. (2) and subtracting Eq. (6) from Eq. (3) yields

(10) {\displaystyle \frac{\partial (u-U)}{\partial t}}={\displaystyle \frac{\mathrm{1}}{D^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \sigma }}\left(Av{\displaystyle \frac{\partial u}{\partial \sigma }}\right)-{\displaystyle \frac{{\tau }_{s\epsilon }-{\tau }_{b\epsilon }}{\rho D}}+M_{\epsilon }-{\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{0}}{M}_{\epsilon }d\sigma },

(11) {\displaystyle \frac{\partial (v-V)}{\partial t}}={\displaystyle \frac{\mathrm{1}}{D^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial \sigma }}\left(Av{\displaystyle \frac{\partial v}{\partial \sigma }}\right)-{\displaystyle \frac{{\tau }_{s\eta }-{\tau }_{b\eta }}{\rho D}}+M_{\eta }-{\displaystyle {\int }_{-\mathrm{1}}^{\mathrm{0}}{M}_{\eta }d\sigma }.

In the above equations, surface gradient terms are eliminated and the relationship between u,v and U,V can be established.

Subtracting Eq. (4) from Eq. (1) yields

(12) {\displaystyle \frac{\partial (h_{\mathrm{2}}D(u-U))}{\partial \epsilon }}+{\displaystyle \frac{\partial (h_{\mathrm{1}}D(v-V))}{\partial \eta }}+h_{\mathrm{1}}{h}_{\mathrm{2}}{\displaystyle \frac{\partial w}{\partial \sigma }}=\mathrm{0.}

If the horizontal velocity is known, then the vertical velocity w can also be obtained through the discretization of Eq. (12).

Boundary conditions

The bottom boundary condition is based on the assumption that the current velocity near the bottom has a logarithmic profile as a function of the distance from the wall. Assume that \overrightarrow{u}(\sigma ) is the horizontal flow rate at a distance of σ from the bottom, which can be expressed by ( Pan et al., 2005)

(13) \overrightarrow{u}(\sigma )={\displaystyle \frac{{\overrightarrow{u}}^{*}}{k}}\ln \{(\mathrm{1}+\sigma )D/z\mathrm{0}\},

where {\overrightarrow{u}}^{*} is the friction velocity, k=0.41 is the Von Karman constant, and z\mathrm{0} is the roughness parameter. The bottom stress is

(14) {\overrightarrow{\tau }}_b=\rho {\overrightarrow{u}}^{*\mathrm{2}}.

By substituting Eq. (13) into Eq. (14), the bottom stress can be expressed as

(15) {\overrightarrow{\tau }}_b=\rho {\left[{\displaystyle \frac{k}{\ln \{(\mathrm{1}+\sigma )D/z\mathrm{0}\}}}\right]}^{\mathrm{2}}\overrightarrow{u}{(\sigma )}^{\mathrm{2}}.

Upon introducing the drag coefficient Cd={\left[{\displaystyle \frac{k}{\ln \{(\mathrm{1}+\sigma )D/z\mathrm{0}\}}}\right]}^{\mathrm{2}} , the bottom stress is given as follows:

(16) ({\tau }_{b\epsilon },{\tau }_{b\eta })=\rho \cdot Cd\sqrt{u_b^{\mathrm{2}}+v_b^{\mathrm{2}}}(u_b,v_b),

where u_b and v_b represent the current velocity nearest to the bottom, satisfying the “logarithmic law” of the velocity profile. Numerically, u_b and v_b are the velocity components of the first grid points nearest to the bottom.

At the free surface (\sigma =\mathrm{0}) , the wind stress components can be expressed as

(17) {\tau }_{sx}={\gamma }_s{\rho }_a{W}_{\epsilon }\sqrt{W_{\epsilon }^{\mathrm{2}}+W_{\eta }^{\mathrm{2}}}\hspace{1em}{\tau }_{sy}={\gamma }_s{\rho }_a{W}_{\eta }\sqrt{W_{\epsilon }^{\mathrm{2}}+W_{\eta }^{\mathrm{2}}},

where W_ε and W_η are the ε and η components of the wind speed vector, ρ_a is the air density; and {\gamma }_s is the wind stress drag coefficient, which is dependent on the wind speed. For the 10-minute-averaged wind speed at 10 m above the free surface, this coefficient can be expressed as ( Garratt, 1977)

(18) {\gamma }_s=\mathrm{6.9}\cdot {\mathrm{10}}^{-\mathrm{4}}+\mathrm{7.5}\cdot {\mathrm{10}}^{-\mathrm{5}}W.

At the open boundaries, the surface water elevations or flow speeds are prescribed as a function of time and space. At the land boundaries, slip conditions are used for the tangential velocity, and a value of zero is prescribed for the normal velocity.

Vertical and horizontal eddy viscosities

The vertical eddy viscosity Av is obtained by introducing appropriate ‘closure’ hypotheses, such as the mixing length hypothesis, the k-ε model, the k-L model, and the eddy viscosity model ( Reid, 1957; Svensson, 1978; Koutitas, 1980; Tsanis, 1989). Performance comparisons of four turbulence closure models implemented via a generic length-scale method are available ( Warner et al., 2005), and there are also several reviews about second-moment turbulence closure models for geophysical boundary layers ( Umlauf and Burchard, 2005; Kantha, 2006). A desirable turbulence sub-model should remain stable and occupy moderate resources during the computation process.

The horizontal eddy viscosity A_H is either simply taken to be a constant value or related to the scales of the motion being resolved, as suggested by Smagorinsky ( Wang and Hutter, 1998)

Numerical scheme

First, we discretize the governing equations into algebraic equations and then find an efficient method of solving them.

Equation discretization

For the discretization of the equations, the staggered C-grid proposed by Arakawa and Lamb is used as shown in Fig. 2. The centers of the cells are numbered with indices i, j, and k, where i=1,…, im in the ε direction, j=1,…, jm in the η direction and k=1,…, km in the vertical direction from the bottom to the free surface. The u velocity is defined at (i+1/2, j, k); The v velocity is defined at (i, j+1/2, k), and the vertical velocity, w, is defined at the node (i, j, k‒1/2). The surface elevation η is defined at the cell center (i, j), and the water depth H is specified at the vertices of each cell to provide a comprehensive representation of the bathymetry.

In accordance with the ADI method, a fractional step method is used to solve the equations in two half steps. In the first half step, the momentum Eqs. (6) and (11) in the η direction and the continuity Eq. (4) are discretized together in the η-σ vertical plane as follows:

(19) {\displaystyle \frac{{\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}}-{\varsigma }_{i,j}^n}{\mathrm{0.5}\Delta t}}{(\Delta h_{\mathrm{1}}\Delta h_{\mathrm{2}})}_{i,j}+{(\Delta h_{\mathrm{2}}{D}^n{U}^n)}_{i+\mathrm{1}/\mathrm{2},j}-{(\Delta h_{\mathrm{2}}{D}^n{U}^n)}_{i-\mathrm{1}/\mathrm{2},j}+{(\Delta h_{\mathrm{1}}{D}^n{V}^{n+\mathrm{1}})}_{i,j+\mathrm{1}/\mathrm{2}} -{(\Delta h_{\mathrm{1}}{D}^n{V}^{n+\mathrm{1}})}_{i,j-\mathrm{1}/\mathrm{2}}=\mathrm{0},

(20) {\displaystyle \frac{V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}-V_{i,j-\mathrm{1}/\mathrm{2}}^n}{\Delta t}}=-g{\displaystyle \frac{{\varsigma }_{i,j}^{n+\mathrm{1}}-{\varsigma }_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}}{{(\Delta h_{\mathrm{2}})}_{i,j-\mathrm{1}/\mathrm{2}}}}+{\displaystyle \frac{{({\tau }_{s\eta })}_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}}{\rho D_{i,j-\mathrm{1}/\mathrm{2}}^n}}-{\displaystyle \frac{C_d\sqrt{{(u_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^n)}^{\mathrm{2}}+{(v_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^n)}^{\mathrm{2}}}}{D_{i,j-\mathrm{1}/\mathrm{2}}^n}}{v}_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}} +{\displaystyle \sum _{k=\mathrm{1}}^{km}{(M_{\eta })}_{i,j-\mathrm{1}/\mathrm{2},k}^n\Delta {\sigma }_k},

(21) {\displaystyle \frac{v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}}-V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}-(v_{i,j-\mathrm{1}/\mathrm{2},k}^n-V_{i,j-\mathrm{1}/\mathrm{2}}^n)}{\Delta t}}={(M_{\eta })}_{i,j-\mathrm{1}/\mathrm{2},k}^n-{\displaystyle \sum _{k=\mathrm{1}}^{km}{(M_{\eta })}_{i,j-\mathrm{1}/\mathrm{2},k}^n\Delta {\sigma }_k}+{\displaystyle \frac{\mathrm{1}}{{(D_{i,j-\mathrm{1}/\mathrm{2}}^n)}^{\mathrm{2}}}}{\displaystyle \frac{A{v}_{i,j-\mathrm{1}/\mathrm{2},k+\mathrm{1}/\mathrm{2}}^n\left({\displaystyle \frac{v_{i,j-\mathrm{1}/\mathrm{2},k+\mathrm{1}}^{n+\mathrm{1}}-v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}}}{\Delta {\sigma }_{k+\mathrm{1}/\mathrm{2}}}}\right)-A{v}_{i,j-\mathrm{1}/\mathrm{2},k-\mathrm{1}/\mathrm{2}}^n\left({\displaystyle \frac{v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}}-v_{i,j-\mathrm{1}/\mathrm{2},k-\mathrm{1}}^{n+\mathrm{1}}}{\Delta {\sigma }_{k-\mathrm{1}/\mathrm{2}}}}\right)}{\Delta {\sigma }_k}}+{\displaystyle \frac{Cd\sqrt{{(u_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^n)}^{\mathrm{2}}+{(v_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^n)}^{\mathrm{2}}}}{D_{i,j-\mathrm{1}/\mathrm{2}}^n}}{v}_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}}-{\displaystyle \frac{{({\tau }_{s\eta })}_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}}{\rho D_{i,j-\mathrm{1}/\mathrm{2}}^n}}.

In the second half step, the process is similar, but the discretization is performed in the ε-σ vertical plane, yielding the following:

(22) {\displaystyle \frac{{\varsigma }_{i,j}^{n+\mathrm{1}}-{\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}}}{\mathrm{0.5}\Delta t}}{(\Delta h_{\mathrm{1}}\Delta h_{\mathrm{2}})}_{i,j}+{(\Delta h_{\mathrm{2}}{D}^n{U}^{n+\mathrm{1}})}_{i+\mathrm{1}/\mathrm{2},j}-{(\Delta h_{\mathrm{2}}{D}^n{U}^{n+\mathrm{1}})}_{i-\mathrm{1}/\mathrm{2},j}+{(\Delta h_{\mathrm{1}}{D}^n{V}^n)}_{i,j+\mathrm{1}/\mathrm{2}} -{(\Delta h_{\mathrm{1}}{D}^n{V}^n)}_{i,j-\mathrm{1}/\mathrm{2}}=\mathrm{0},

(23) {\displaystyle \frac{U_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}}-U_{i-\mathrm{1}/\mathrm{2},j}^n}{\Delta t}}=-g{\displaystyle \frac{{\varsigma }_{i,j}^{n+\mathrm{1}}-{\varsigma }_{i-\mathrm{1},j}^{n+\mathrm{1}}}{{(\Delta h_{\mathrm{1}})}_{i-\mathrm{1}/\mathrm{2},j}}}+{\displaystyle \frac{{({\tau }_{s\epsilon })}_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}}}{\rho D_{i-\mathrm{1}/\mathrm{2},j}^n}}-{\displaystyle \frac{C_d\sqrt{{(u_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^n)}^{\mathrm{2}}+{(v_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^n)}^{\mathrm{2}}}}{D_{i-\mathrm{1}/\mathrm{2},j}}}{u}_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^{n+\mathrm{1}} +{\displaystyle \sum _{k=\mathrm{1}}^{km}{(M_{\epsilon })}_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}/\mathrm{2}}\mathrm{\Delta }{\sigma }_k},

(24) {\displaystyle \frac{u_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}}-U_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}}-(u_{i-\mathrm{1}/\mathrm{2},j,k}^n-U_{i-\mathrm{1}/\mathrm{2},j}^n)}{\mathrm{\Delta }t}}={(M_{\epsilon })}_{i-\mathrm{1}/\mathrm{2},j,k}^n-{\displaystyle \sum _{k=\mathrm{1}}^{km}{(M_{\epsilon })}_{i-\mathrm{1}/\mathrm{2},j,k}^n\mathrm{\Delta }{\sigma }_k}+{\displaystyle \frac{\mathrm{1}}{{(D_{i-\mathrm{1}/\mathrm{2},j}^n)}^{\mathrm{2}}}}{\displaystyle \frac{A{v}_{i-\mathrm{1}/\mathrm{2},j,k+\mathrm{1}/\mathrm{2}}^n\left({\displaystyle \frac{u_{i-\mathrm{1}/\mathrm{2},j,k+\mathrm{1}}^{n+\mathrm{1}}-u_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}}}{\mathrm{\Delta }{\sigma }_{k+\mathrm{1}/\mathrm{2}}}}\right)-A{v}_{i-\mathrm{1}/\mathrm{2},j,k-\mathrm{1}/\mathrm{2}}^n\left({\displaystyle \frac{u_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}}-u_{i-\mathrm{1}/\mathrm{2},j,k-\mathrm{1}}^{n+\mathrm{1}}}{\mathrm{\Delta }{\sigma }_{k-\mathrm{1}/\mathrm{2}}}}\right)}{\mathrm{\Delta }{\sigma }_k}}+{\displaystyle \frac{Cd\sqrt{{(u_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^n)}^{\mathrm{2}}+{(v_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^n)}^{\mathrm{2}}}}{D_{i-\mathrm{1}/\mathrm{2},j}^n}}{u}_{i-\mathrm{1}/\mathrm{2},j,\mathrm{1}}^{n+\mathrm{1}}-{\displaystyle \frac{{({\tau }_{s\epsilon })}_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}}}{\rho D_{i-\mathrm{1}/\mathrm{2},j}^n}},

where {\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}} denotes the intermediate water elevation at time t_n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t and \mathrm{\Delta }{\sigma }_{k\pm \mathrm{1}/\mathrm{2}}=\mathrm{0.5}\times (\mathrm{\Delta }{\sigma }_k+\mathrm{\Delta }{\sigma }_{k\pm \mathrm{1}}) is the vertical increment between levels k and k±1. Here, we define \Delta h_{\mathrm{1}}=h_{\mathrm{1}}\Delta \epsilon and \Delta h_{\mathrm{2}}=h_{\mathrm{2}}\Delta \eta , which are the horizontal grid scales in both directions.

Computational method

Taking the first half step as an example, Eq. (21) can be rewritten as

(25) A{\mathrm{1}}_{k+\mathrm{1}}{v}_{i,j-\mathrm{1}/\mathrm{2},k+\mathrm{1}}^{n+\mathrm{1}}+A{\mathrm{2}}_k{v}_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}}+A{\mathrm{3}}_{k-\mathrm{1}}{v}_{i,j-\mathrm{1}/\mathrm{2},k-\mathrm{1}}^{n+\mathrm{1}}+A\mathrm{4}{v}_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}}+A\mathrm{5}{V}_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}=A{\mathrm{6}}_k,

where A{\mathrm{1}}_{k+\mathrm{1}} , A{\mathrm{2}}_k , A{\mathrm{3}}_{k-\mathrm{1}} , A\mathrm{4} , A\mathrm{5} and A{\mathrm{6}}_k are coefficients obtained at the last time level.

An auxiliary equation at level k is introduced as follows:

(26) v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}}=P_k{V}_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}+Q_k{v}_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}}+R_k.

It is obvious that P_{\mathrm{1}}=\mathrm{0} , Q_{\mathrm{1}}=\mathrm{1} and R_{\mathrm{1}}=\mathrm{0} .

In addition, at level k=1, P_{\mathrm{2}}=-A\mathrm{5}/A{\mathrm{1}}_{\mathrm{2}} , Q_{\mathrm{2}}=-(A\mathrm{4}+A{\mathrm{2}}_{\mathrm{1}})/A{\mathrm{1}}_{\mathrm{2}} and R_{\mathrm{2}}=A{\mathrm{6}}_{\mathrm{1}}/A{\mathrm{1}}_{\mathrm{2}} are obtained from Eq. (25).

Substituting Eq. (26) at levels k=m‒1 and m (m>1) into Eq. (25) yields

(27) v_{i,j-\mathrm{1}/\mathrm{2},m+\mathrm{1}}^{n+\mathrm{1}}=-{\displaystyle \frac{A\mathrm{5}+A{\mathrm{2}}_m{P}_m+A{\mathrm{3}}_{m-\mathrm{1}}{P}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}}{V}_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}-{\displaystyle \frac{A\mathrm{4}+A{\mathrm{2}}_m{Q}_m+A{\mathrm{3}}_{m-\mathrm{1}}{Q}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}}{v}_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}}+{\displaystyle \frac{A{\mathrm{6}}_m-A{\mathrm{2}}_m{R}_m-A{\mathrm{3}}_{m-\mathrm{1}}{R}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}}.

A comparison between Eqs. (26) and (27) reveals that P_{m+\mathrm{1}} , Q_{m+\mathrm{1}} and R_{m+\mathrm{1}} (m>1) can be expressed as

(28a) P_{m+\mathrm{1}}=-{\displaystyle \frac{A\mathrm{5}+A{\mathrm{2}}_m{P}_m+A{\mathrm{3}}_{m-\mathrm{1}}{P}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}},

(28b) Q_{m+\mathrm{1}}=-{\displaystyle \frac{A\mathrm{4}+A{\mathrm{2}}_m{Q}_m+A{\mathrm{3}}_{m-\mathrm{1}}{Q}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}},

(28c) R_{m+\mathrm{1}}={\displaystyle \frac{A{\mathrm{6}}_m-A{\mathrm{2}}_m{R}_m-A{\mathrm{3}}_{m-\mathrm{1}}{R}_{m-\mathrm{1}}}{A{\mathrm{1}}_{m+\mathrm{1}}}}.

Finally, given a total number of km layers in the vertical direction, multiplying Eq. (27) by \mathrm{\Delta }{\sigma }_k , which is the relative thickness of depth layer k, and then summing from k=1 to K yields

(29) v_{i,j-\mathrm{1}/\mathrm{2},\mathrm{1}}^{n+\mathrm{1}}=\left[\left(\mathrm{1}-{\displaystyle \sum _{k=\mathrm{1}}^K{P}_k\mathrm{\Delta }{\sigma }_k}\right)V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}-{\displaystyle \sum _{k=\mathrm{1}}^K{R}_k\mathrm{\Delta }{\sigma }_k}\right]/{\displaystyle \sum _{k=\mathrm{1}}^K{Q}_k\mathrm{\Delta }{\sigma }_k}.

Then, upon substituting Eq. (29) into Eq. (20), the surface elevation {\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}} and the depth-averaged velocity component V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}} can be solved for implicitly using the ADI method; thus, the velocity profiles can be obtained simultaneously from Eqs. (29) and (26). For the next half step, the process is the same as that for the first, except with v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}} and V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}} replaced by u_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}} and U_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}} .

Summary of the scheme

In summary, each time step is divided into two half steps. In the first half step, the coefficients P_i, Q_i, and R_i (i=1, 2, …, km) are obtained from Eq. (21). Then, Eq. (29) should be substituted into Eq. (20). Afterward, Eqs. (19) and (20) are solved, combined with the boundary conditions for {\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}} and V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}} . Using Eqs. (29) and (26), v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}} can be obtained. Finally, the vertical velocities w_{i,j,k}^{n+\mathrm{1}} are obtained through the discretization of Eq. (12). The process in the second half time step is much like that in the first, except that the direction changes from η to ε and {\varsigma }_{i,j}^n , v_{i,j-\mathrm{1}/\mathrm{2},k}^{n+\mathrm{1}} and V_{i,j-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}} are replaced with {\varsigma }_{i,j}^{n+\mathrm{1}/\mathrm{2}} , u_{i-\mathrm{1}/\mathrm{2},j,k}^{n+\mathrm{1}} and U_{i-\mathrm{1}/\mathrm{2},j}^{n+\mathrm{1}} .

Wetting and drying (WAD) scheme

Regarding the tidal model, an important phenomenon that should be addressed using a specialized numerical method is wetting and drying (WAD). This phenomenon typically occurs in low-lying coastal zones and also in embayments and inlets. This process should be reflected in any tidal model, such as the well-known POM ( Oey, 2005; Xue and Du, 2010). At present, several WAD evaluations have been performed ( Medeiros and Hagen, 2013). Each WAD algorithm is closely related to the numerical model it serves, and thus, they are difficult to generalize.

Overall, a robust WAD scheme not only guarantees stability and mass balance but also reflects the real physical process precisely. In this paper, the improved line boundary method is used ( Song et al., 1999). In this method, the side of the cell rather than the entire cell is used to judge the wetting and drying status. If the water depth at one side of the cell is lower than a certain threshold value, then that side is treated as a closed boundary, called a line boundary. In one grid cell, the minimum number of line boundaries is zero and the maximum is five; different combinations of line boundaries can represent different boundary forms. Thus, this method can simulate the wetting and drying process more accurately than can the conventional method.

Validation of the model

A crucial step in numerical modeling is the verification of the model. Three test cases are presented in this paper. First, the model is validated against the analytical solution developed by Lynch (1985), which is often used in the literature to test numerical solutions ( Toro and Gomez, 1998).

The geometry and topography of the test case are the same as those of the case used by Toro and Gomez, shown in Fig. 3. It is an annular bay with a depth following a parabolic distribution h=h_{\mathrm{0}}{r}^{\mathrm{2}} ( h_{\mathrm{0}} is a constant) from the shallower boundary at r1 to the deeper one at r2. A semi-diurnal tidal wave of variable amplitude is prescribed along the boundary at r\mathrm{2} . The bottom stress is assumed to vary linearly with the velocity, and the vertical eddy viscosity coefficient is assumed to be constant. The parameters used in the test are as follows: r1=6096.0 m, r2=152400.0 m, h₀=3.048r1⁻², T=12 hours, {\zeta }_{r\mathrm{2}}=\cos (\mathrm{2}\theta ) for 0≤θ≤ \pi /2.

Considering that the analytical solution is obtained from linearized mass and momentum conservation equations without convective and horizontal diffusive terms, we leave out these terms in the numerical analysis for comparison with the analytical result.

The model was run for several periods, until the establishment of steady periodic circulation. To test the stability and sensitivity of the model with different Courant numbers and vertical resolutions, several time steps (i.e., 400 s, 200 s, 100 s, 50 s) and vertical levels (i.e., 5, 8, 11) were chosen for comparison with each other. It was found that even at large Courant numbers, when a relatively large time step and high vertical resolution were used, the model was stable and there was no obvious effect on the current pattern and tide level. Figure 4 shows the free surface elevation contours and the velocity field at the surface in various phases of a single period.

The free surface lines along sections 1-1, 2-2, and 3-3 are shown in Fig. 5, and the vertical profiles of the x and y velocity components are shown in Fig. 6. The agreement between the analytical and numerical solutions is very good, with the exception of a few small differences, possibly due to the inversion of the velocity from curvilinear coordinates to Cartesian coordinates using a relatively coarse grid and the numerical diffusion caused by the scheme itself. Overall, the result is satisfactory.

The second test case is a simulation of a wind-driven current in an idealized rectangular basin. The vertical profile of the velocity is available from Officer (1976) and is written as follows:

(30) {\displaystyle \frac{\partial \zeta }{\partial x}}={\displaystyle \frac{{\displaystyle \frac{{\tau }_s}{\rho }}\left({\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{2}Av}}+{\displaystyle \frac{h}{{\gamma }_b}}\right)}{\left({\displaystyle \frac{h^{\mathrm{3}}}{\mathrm{3}Av}}+{\displaystyle \frac{h^{\mathrm{2}}}{{\gamma }_b}}\right)}},

(31) u(z)=g{\displaystyle \frac{\partial \zeta }{\partial x}}\left({\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}Av}}-{\displaystyle \frac{h}{{\gamma }_b}}-{\displaystyle \frac{h^{\mathrm{2}}}{\mathrm{2}Av}}\right)+{\displaystyle \frac{{\tau }_s}{\rho }}\left({\displaystyle \frac{z}{Av}}+{\displaystyle \frac{\mathrm{1}}{{\gamma }_b}}+{\displaystyle \frac{h}{Av}}\right),

where {\gamma }_b is the bed friction coefficient and is assumed to be 0.001 m/s.

In the calculation, the rectangular closed basin is 3,400 m long, 1,400 m wide, and 10 m in depth. The wind blows longitudinally to the basin, and a constant vertical turbulent eddy coefficient of A_v=\mathrm{0.01} is assumed. The other parameters are chosen to be the same as in the test case considered by Drago and Iovenitti ( 2000).

A comparison between the numerical and analytical vertical current profiles in the middle of the basin is shown in Fig. 7 for different numbers of layers at a wind speed of 22.8 m/s. The numerical results can be improved by using a higher vertical resolution. It can be observed that the action of the wind on the water’s surface causes a drift along the direction of the wind and a bottom return flow in the opposite direction. This indicates that the model can reproduce the circulation for a geophysical wind-driven current.

The third test is to investigate the model’s ability to reproduce the Ekman spiral. As we know, the Coriolis effect plays a major role in coastal dynamics, and its correct representation is a critical requirement for large-scale tidal models. Ekman dynamics describes the expected behavior for a steady wind blowing over an infinitely deep and wide ocean with a constant density, assuming a balance among the wind stress, vertical eddy viscosity, and the Coriolis force. Under the assumption that the wind is blowing longitudinally in the easterly direction, the analytical solution is given by

(32) \begin{array}{l}u=V_{\mathrm{0}}{e}^{-(\pi /D_{\mathrm{0}})z}\sin \left({\displaystyle \frac{\pi }{\mathrm{4}}}-{\displaystyle \frac{\pi }{D_{\mathrm{0}}}}z\right)\hfill \\ v=V_{\mathrm{0}}{e}^{-(\pi /D_{\mathrm{0}})z}\cos \left({\displaystyle \frac{\pi }{\mathrm{4}}}-{\displaystyle \frac{\pi }{D_{\mathrm{0}}}}z\right)\hfill \end{array}.

where D_{\mathrm{0}} is the Ekman depth:

(33) D_{\mathrm{0}}=\pi \sqrt{{\displaystyle \frac{Av}{\omega \sin \mathrm{\Phi }}}}.

This parameter represents the height of the water layer influenced by the wind, and V_{\mathrm{0}} is the surface current speed:

(34) V_{\mathrm{0}}={\displaystyle \frac{{\tau }_s}{\rho }}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}Av\omega \sin \phi }}},

where \omega is the angular velocity of the Earth and \phi is the latitude.

The analytical solution shows the currents at an angle to the direction of the wind: π/4 to the right of the wind direction at the surface in the northern hemisphere, and rotating with depth in a spiral pattern known as the Ekman spiral. The computational domain, a flat layer of 100 km×100 km with a frictionless bottom (to approximate an infinite depth) at a latitude of approximately 45°, is discretized using 2 km×2 km meshes and 15 σ layers with thicknesses of 0.01 (10 layers), 0.1 (l layer), and 0.2 (4 layers) from the surface to the bottom. The elevation is fixed at the mean sea level (MSL) on the four edges of the square domain. The time step is set to 600 s. The non-linear advection terms are turned off to facilitate comparison with the analytical solution. The other parameters used for the simulation are as follows: Av=0.01 m²/s and τ_s=0.176 N/m² (when the wind speed is 10 m/s to the east).

The solution oscillates around its steady state with a velocity amplitude of approximately 0.05 m/s for a long time. To avoid excessive CPU time consumption, we average the solution at the center of the domain over more than 1 day to approximate a quasi-steady state. Figure 8 shows a comparison between the analytical and computed solutions for the Ekman spiral.

To evaluate the performance of the model and demonstrate its predictive capacity in a practical case, we consider its application to the Jiangsu offshore sea area, located to the east of Jiangsu province in China in a near-shore region of the East China Sea, which is morphologically complex, especially near the radiation shoal that has formed under specific hydrodynamic conditions for more than a thousand years and is considered to have a significant effect on the tidal wave system in the East China Sea. As shown in Fig. 9, the simulation domain extends north to the southern coastline of Shandong Peninsula and south to the Yangtze River Estuary, spanning approximate dimensions of 630 km in length and 106 km in width. Its northern and western boundaries are coastlines, and its eastern and southern boundaries are open ocean and approximately perpendicular.

The horizontal grid employs a curvilinear orthogonal system with a variable resolution, ranging from 1−0.5 km near the open sea in the east to 0.2−0.5 km at the Yangtze River Estuary. The vertical σ grid has 7 levels, with higher resolutions toward the bottom.

Lateral open boundary conditions for the tide level are provided by the East China Sea Model (Song et al., 1999). The normal vector of the velocity is defined to be zero at land boundaries. The upstream water level of the Yangtze River is supplied by hydrologic stations.

The parameters used in the model are as follows: a time interval of 100 s, Ah=100 m²/s, and z0=0.01 cm.

The numerical results are compared with the water levels and velocities measured at select measurement sites distributed at monitoring stations a, b, c, and d (see Fig. 10 and Fig. 11). Stations a, b, c, and d were chosen for water-level comparisons, and Stations b and c, including 7 and 5 sites, respectively, were chosen for current velocity comparisons. Here, we show only some of the sites for the current velocity comparisons over more than 1 tidal cycle. The predicted tide level is consistent with the measured one in terms of both phasing and amplification. Note, however, that the largest discrepancies appear mainly at wave crests and troughs. In addition to inadequate model physics, these discrepancies are, in part, due to the inexact knowledge of the open boundary conditions resulting from the linear interpolation from the coarse grids of the East China Sea Model.

The red points represent the measurement sites used to validate the numerical simulation results and the blue area represents the simulation domain. The right panel shows the sand ridge of the radiation shoal.

After model verification, we can determine the factor for the formation of the radiation shoal in terms of the hydrodynamic angle. Figure 12 shows the surface velocity field at times of peak flood and peak ebb. It is observed that the tide current is characteristic of periodic to-and-fro movement. The current between the Yangtze River Estuary and the radiation shoal rotates clockwise, and the current in the north near Shandong Peninsula rotates anticlockwise. When the current floods, it flows radiatively into the offshore region, and the tidal bores converge on the radiation shoal in a ring-like shape. When the current ebbs, it flows radiatively out of the offshore region, and the tidal bores diverge away from the radiation shoal. This unique flow pattern is caused by the joint effect of tide waves traveling from the East China Open Sea and the particular coastline and continent shelf geometry of the Shandong Peninsula and Jiangsu Province. The current structure plays a major role in the evolution of the radiation shoal and the morphology of the offshore water. The water depth near the radiation shoal is shallow; therefore, large areas of sand flats are exposed before the tidal bores pass over them. The WAD schemes applied in the model can effectively address this problem.

As another example, we demonstrate the application of the model to the tidal dynamics system of the Bohai Sea and the Yellow Sea. As shown in Fig. 13, the computational domain is discretized using a longitude-latitude grid, which is a type of orthogonal curvilinear grid. The spatial resolution is set to 2 minutes (approximately 0.333 degree) along each direction. The southern and south-eastern boundaries of the domain are specified as open boundaries, forced by 8 major constituents, 4 semi-diurnal and 4 diurnal, which are typically written as M₂, S₂, N₂, K₂, K₁, O₁, P₁, and Q₁. We ran the model at a time step of 450 s for 33 mode-days, and the time series of surface elevations at each grid point was then fitted using the least-squares harmonic analysis method.

To test the accuracy of our model, 10 tidal stations along the coastline were selected, of which the last 3 stations are located on the western shoreline of the Korean peninsula (see the left-subfigure in Fig. 13). The harmonic constants obtained from these stations and from our model are listed in Table 1, where the first 2 columns give the name and location of each station and the last 4 columns show a comparison of the harmonic constants for constituents M₂, S₂, K₁, and O_1. From the overall root-mean-square (RMS) differences, it can be observed that the model offers greater accuracy in amplitude than in phase lag. Moreover, comparisons between the harmonic constants obtained from OPEX/Poseidon altimetry ( Fang, et al., 2004) and from the numerical model are given in Table 2. The RMS differences in amplitude for M₂, S₂, K₁, and O₁ are less than 8 cm, and the differences in phase lag are not more than 10°, indicating that the difference in accuracy between the two models is small and may be tolerable.

The resulting co-tidal charts for each of these four constituents are shown in Fig. 14 and Fig. 15. It can be observed from Fig. 14 that the M₂ and S₂ tides have three amphidromic points: one in the Bohai Sea and two in the Yellow Sea. The amphidromic point in the Bohai Sea is close to the mouth of the Yellow River. The two amphidromic points in the Yellow Sea are located to the northeast of Chengshantou and the southeast of Qingdao. The pattern of the S₂ constituent is much like that of M₂, differing only slightly in the details. Figure 15 shows that the K₁ and O₁ constituents have two amphidromic points: one point in the Bohai Strait and the other in the Southern Yellow Sea. Similarly, the pattern of the O₁ constituent is much like that of K_1. These results agree with those obtained by other authors ( Kang et al., 1998; Fang et al., 2004).

Conclusions

A 3D tidal model has been established that can be applied to estuarine, coastal, and near-shore areas with rather complex geometry and bathymetry. A terrain-following σ coordinate system is applied vertically to represent the bathymetry as accurately as possible. An orthogonal curvilinear coordinate system is incorporated into the model to better fit the coastline. Another characteristic of the model is the numerical scheme, which is based on the conventional ADI method; it inherits certain advantages of this method, such as high efficiency and stability, and can be extended to 3D problems with little computational effort.

The model was tested against several cases, including analytical solutions and measured data. In all cases, the model yields very satisfactory results, thus confirming the accuracy and efficiency of the model for engineering applications. In addition, the model was applied to the Jiangsu offshore sea area, a near-shore region of the Yellow Sea to the east of China. The model was also used to predict the patterns of the tidal wave system in the Bohai Sea and Yellow Sea. These applications and simulations will be helpful for studying the dynamics of a tidal wave system and its interactions with geomorphology.

Further developments of the model will account for the baroclinic effect, achieve coupling with a turbulence sub-model, and incorporate a pollutant or sediment transport module. In addition, the numerical model should be applied to more practical cases for calibration and improvement.

References

[1] Bell M J (1999). Vortex stretching and bottom torque in the Bryan-Cox ocean circulation model. J Geophys Res, 104(C10): 23545–23563

[2] Blumberg A F, Herring H J (1987). Circulation modelling using orthogonal curvilinear coordinates. In: Nihoul J C J, Jamart B M, eds. Three-dimensional Models of Marine and Estuarine Dynamics. Elsevier Oceanography Series, (45): 221–243

[3] Blumberg A F, Mellor G L (1987). A description of a three-dimensional coastal ocean circulation model. In: Heaps N ed. Three-dimensional Coastal Ocean Models. Coastal and Estuarine Sciences, 4: 1–16

[4] Casulli V, Stelling G S (1998). Numerical simulation of 3D quasi-hydrostatic, free-surface flows. J Hydraul Eng, 124(7): 678–686

[5] Chen Q, Tan K, Zhu C, Li R (2009). Development and application of a two-dimensional water quality model for the Daqinghe River Mouth of the Dianchi Lake. J Environ Sci (China), 21(3): 313–318

[6] Drago M, Iovenitti L (2000). σ-Coordinates hydrodynamic numerical model for coastal and ocean three-dimensional circulation. Ocean Eng, 27(10): 1065–1085

[7] Fang G H, Wang Y G, Wei Z X, Choi B H, Wang X Y, Wang J (2004). Empirical cotidal charts of the Bohai, Yellow and East China Seas from 10 years of TOPEX/Poseidon altimetry. J Geophys Res, 109(C11): C11006

[8] Garratt J R (1977). Review of drag coefficients over oceans and continents. Mon Weather Rev, 105(7): 915–929

[9] Gerdes R (1993). A primitive equation ocean circulation model using a general vertical coordinate transformation. Part 1: description and testing of the model. J Geophys Res, 98(C8): 14683–14701

[10] Hu D C, Fan B L, Wang G Q, Zhang H W (2011a). A semi-implicit 3-D numerical model using sigma-coordinate for non-hydrostatic pressure free-surface flows. J Hydrodynam, 23(2): 212–223

[11] Hu K L, Ding P X, Wang Z B, Yang S L (2009). 2D/3D hydrodynamic and sediment transport model for the Yangtze Estuary, China. J Mar Syst, 77(1−2): 114–136

[12] Hu K M, Pang Y, Wang H, Wang X M, Wu X W, Bao K, Liu Q (2011b<?Pub Caret?>). Simulation study on water quality based on sediment release flume experiment in Lake Taihu, China. Ecol Eng, 37(4): 607–615

[13] Kang S K, Lee S R, Lie H J (1998). Fine grid tidal modeling of the Yellow and East China Seas. Cont Shelf Res, 18(7): 739–772

[14] Kantha L (2006). Comments on ‘‘Second-order turbulence closure models for geophysical boundary layers: a review of recent work’’. Cont Shelf Res, 26(6): 819–822

[15] Keilegavlen E, Berntsen J (2009). Non-hydrostatic pressure in σ coordinate ocean models. Ocean Model, 28(4): 240–249

[16] Koutitas C (1980). Modelling three-dimensional wind-induced flows. J Hydraul Div, 106(11): 1843–1865

[17] Lazure P, Dumas F(2008). An external−internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Adv Water Resour, 31(2): 233–250

[18] Leendertse J J (1970). Water quality model for well-mixed estuary and coastal seas, Vol.1 Principle of computation. The Rand Corporation, Report No. RM-6230-RC, Santa Monica, California

[19] Ly L N, Luong P (1999). Numerical multi-block grids in coastal ocean circulation modeling. Appl Math Model, 23(11): 865–879

[20] Mao J, Chen Q, Chen Y (2008). Three-dimensional eutrophication model and application to Taihu Lake, China. J Environ Sci (China), 20(3): 278–284

[21] Medeiros S C, Hagen S C (2013). Review of wetting and drying algorithms for numerical tidal flow models. Int J Numer Methods Fluids, 71(4): 473–487

[22] Oey L (2005). A wetting and drying scheme for POM. Ocean Model, 9(2): 133–150

[23] Officer C B (1976). Physical Oceanography of Estuaries (and Associated Coastal Waters). New York: John Wiley and Sons, 465

[24] Pan J Y, Wang D W, Hwang P A (2005). A study of wave effects on wind stress over the ocean in a fetch-limited case. J Geophys Res, 110(C2): C02020

[25] Phillips N A (1957). A coordinate system having some special advantages for numerical forecasting. J Meteorol, 14(2): 184–185

[26] Reid R O (1957). Modification of the quadratic bottom-stress law of turbulent channel flow in the presence of surface wind-stress. US Department of the Army, Tech.Mem., 93, Beach Erosion Board, Washington, D.C.

[27] Shi F Y, Ding P X, Kong Y Z (1999). Numerical tidal current modeling using fine boundary-fitted grids for tidal flats. China Ocean Engineering, 30(2): 115–124 (in Chinese)

[28] Song Z Y, Xue H C, Yan Y X, Mao L H, Xu F M (1999). Quasi-3D numerical simulation of tidal hydrodynamic field. China Ocean Engineering, 13(3): 265–276 (in Chinese)

[29] Stansby P K, Zhou J G (1998). Shallow water flow solver with non-hydrostatic pressure: 2D vertical plane problems. Int J Numer Methods Fluids, 28(3): 541–563

[30] Svensson U (1978). Mathematical model of the seasonal thermocline. Report No.1002, Department of Water Resources Engineering, University of Lund, Sweden

[31] Toro B, Gomez G (1998). A semi-implicit 3D numerical model for free surface flows. Parallel Session. Numerical Methods for Current and Wave Calculation, 56: 1–26

[32] Tsanis I (1989). Simulation of wind-induced water currents. J Hydraul Eng, 115(8): 1113–1134

[33] Umlauf L, Burchard H (2005). Second-order turbulence closure models for geophysical boundary layers. A review of recent work. Cont Shelf Res, 25(7−8): 795–827

[34] Wang Y, Hutter K (1998). A semi-implicit semispectral primitive equation model for lake circulation dynamics and its stability performance. J Comput Phys, 139(1): 209–241

[35] Warner J C, Sherwood C R, Arango H G, Signell R P(2005). Performance of four turbulence closure models implemented using a generic length scale method. Ocean Model, 8(1−2): 81–113

[36] Xie M X, Zhuang W (2010). Numerical study on the three-dimensional characteristics of the tidal current around harbor entrance. J Hydrodynam, 22(6): 847–855

[37] Xue H J, Du Y (2010). Implementation of a wetting-and-drying model in simulating the Kennebec−Androscoggin plume and the circulation in Casco Bay. Ocean Dyn, 60(2): 341–357

[38] Yang L W, Ozer J (1997). An ADI technique for solving three-dimensional coastal circulation. Acta Oceanol Sin, 16(3): 281–302

[39] Young C C, Wu C H, Kuo J T, Liu W C (2007). A higher-order σ coordinate non-hydrostatic model for nonlinear surface waves. Ocean Eng, 34(10): 1357–1370

[40] Zhang J X, Liu H, Xue L P (2006). A vertical 2D mathematical model for hydrodynamic flows with free surface in σ coordinate. J Hydrodynam, 18(1): 82–90

[41] Zhang X B, Hu D C, Wang M (2010). A 2-D hydrodynamic model for the river, lake and network system in the JingJiang Reach on the unstructured quadrangles. J Hydrodynam, 22(3): 419–429

[42] Zhao L, Zhang X, Liu Y, He B, Zhu X, Zou R, Zhu Y (2012). Three-dimensional hydrodynamic and water quality model for TMDL development of Lake Fuxian, China. J Environ Sci (China), 24(8): 1355–1363

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg