Simulation of Ocean Circulation of Dongsha Water Using Non-Hydrostatic Shallow-Water Model

Abstract

A two-dimensional non-hydrostatic shallow-water model for weakly dispersive waves is developed using the least-squares finite-element method. The model is based on the depth-averaged, nonlinear and non-hydrostatic shallow-water equations. The non-hydrostatic shallow-water equations are solved with the semi-implicit (predictor-corrector) method and least-squares finite-element method. In the predictor step, hydrostatic pressure at the previous step is used as an initial guess and an intermediate velocity field is calculated. In the corrector step, a Poisson equation for the non-hydrostatic pressure is solved and the final velocity and free-surface elevation is corrected for the new time step. The non-hydrostatic shallow-water model is verified and applied to both wave and flow driven fluid flows, including solitary wave propagation in a channel, progressive sinusoidal waves propagation over a submerged bar, von Karmann vortex street, and ocean circulations of Dongsha Atolls. It is found hydrostatic shallow-water model is efficient and accurate for shallow water flows. Non-hydrostatic shallow-water model requires 1.5 to 3.0 more cpu time than hydrostatic shallow-water model for the same simulation. Model simulations reveal that non-hydrostatic pressure gradients could affect the velocity field and free-surface significantly in case where nonlinearity and dispersion are important during the course of wave propagation.

1. Introduction

An accurate prediction of multi-scale ocean waves and flows is important in physical oceanography, marine hydrodynamics, as well as ocean and coastal engineering. The hydrostatic shallow-water models (SWMs) have been widely used in studying flows in rivers, lakes, estuaries, and oceans [1,2,3]. The hydrostatic assumption is valid in many cases and hydrostatic SWMs have been successfully used in numerous engineering applications. However, there are many cases this assumption is questionable. Due to the hydrostatic and non-dispersive assumption of the hydrostatic SWMs, their application to short period waves, abrupt changes in bed topographies, and stratification due to strong density gradients become limited and result in inaccurate predictions [4,5,6,7,8].

Generally, two groups of numerical methods for relaxing the limitations of the hydrostatic and non-dispersive assumption and simulating deep water waves have been developed. One is the Boussinesq-type model attributed to the pioneer work of [9,10], and another approach is proposed by using the Reynolds-Averaged Navier–Stokes (RANS) equations with the non-hydrostatic pressure [11,12,13,14,15,16,17]. The classical Boussinesq-type equation [9] mainly described vertical velocity profile with a polynomial for weakly nonlinear dispersive waves. Subsequently, the constrains have been further released using optimal expansion levels, higher-degree polynomials, higher-order derivatives, and multiple layers [10,18,19,20].

Many non-hydrostatic models that solve non-overturning free-surface flows using single-valued height method have been developed. Conventionally, a large number of vertical layers (e.g., 10–40) with finer resolution near the free surface are utilized to resolve dispersive and non-linear features of water waves [21,22,23,24,25]. Note that non-hydrostatic pressure is taken into account for most part of the water column but a hydrostatic pressure relation is still adopted at the top-layer cell. Later, efficient non-hydrostatic models [26,27,28,29,30,31,32] rigorously treat the top-layer non-hydrostatic pressure effect. The required vertical layers are dramatically reduced to a relatively small number (2–5 layers, an order of magnitude smaller than those required in traditional models), effectively decreasing the computational expenses.

For example, ref. [16] proposed an edged-based Keller-box scheme to exactly assign zero surface pressure without any approximation. Further applications of this model can be found in [17,32]. For a staggered grid system, ref. [27] utilized an integration method to tackle the top-layer hydrostatic pressure assumption. Besides, several interpolation or extrapolation approaches [28,29,30,31] have also been implemented in non-hydrostatic models. Generally, these models with two or three layers can resolve wave dispersion up to a dimensionless relative water depth kh ≈ 1 or 2.5, respectively, where k and h are the wave number and water depth. Continuous research effort is to establish a theoretical framework for the interpretation of wave dispersion/nonlinearity [32,33] and to predict a broader range of surface wave motions with improved numerical accuracy and efficiency, e.g., [29,34,35,36,37,38]. In addition to shallow-water flows, further application of efficient non-hydrostatic models in predicting near-shore wave transformation and weakly dispersive tsunami propagation [20,21,22,23,24,25,26,27,28,29,30,31,32], therefore, becomes promising.

To date, the non-hydrostatic models have been regarded as the most general simulation tool for the coastal research, especially with the aid of advanced computer capacity today [33,35]. Non-hydrostatic models are often developed based on their hydrostatic predecessors [16,26,28]. In this study, a non-hydrostatic shallow-water model is developed and applied to simulate ocean circulations of Dongsha waters. The model is based on the hydrostatic predecessor [39,40] and employs the semi-implicit method as well as least-squares finite-element formulation. Compared with Galerkin and other numerical methods, least-squares finite-element method enjoys many programming and computational advantages, such as the use of single approximating space for all variables, and the choice of approximating spaces not being subjected to the Ladyzhenskaya–Babuska–Brezzi (LBB) condition. Therefore, mixed formulation, collocation grid arrangement, and convection scheme are unnecessary [41,42]. More importantly, the resulting system of equations of LSFEM is symmetric, positive-definite (SPD), and can be used with an efficient element-to-element conjugate gradient method [43].

Several benchmark cases with available analytical solutions, experimental data, and numerical results from other models are used to validate the model. Finally, the model is applied to simulate ocean circulations of Dongsha waters. The paper is structured as follows. In Section 2, the non-hydrostatic shallow-water equations and least-squares finite-element formulation are introduced. Model verification and applications are presented in Section 3. Based on the computation results, discussions are made in Section 4 and conclusions are drawn in Section 5.

2. Materials and Methods

Depth-integrated two-dimensional non-hydrostatic shallow-water equations and numerical methods are briefly introduced in this section.

2.1. Shallow-Water Equations

The two-dimensional depth-averaged non-hydrostatic shallow-water equations are derived from the depth-averaged Navier–Stokes equations and read.

$$\frac{\partial \eta }{\partial t}+\frac{\partial \left(HU\right)}{\partial x}+\frac{\partial \left(HV\right)}{\partial y}=0$$

(1)

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}=-g\frac{\partial \eta }{\partial x}-\frac{1}{2\rho H}\left(H\frac{\partial q}{\partial x}+q\left(\frac{\partial \eta }{\partial x}+\frac{\partial h}{\partial x}\right)\right)+S_x$$

(2)

$$\frac{\partial U}{\partial t}+U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}=-g\frac{\partial \eta }{\partial y}-\frac{1}{2\rho H}\left(H\frac{\partial q}{\partial y}+q\left(\frac{\partial \eta }{\partial y}+\frac{\partial h}{\partial y}\right)\right)+S_y$$

(3)

$$\frac{DW}{Dt}=\frac{q}{\rho H}$$

(4)

where H(x, y, t) = η(x, y, t) − h(x, y) is the total water depth, h(x, y) is the bottom elevation, η(x, y, t) is the surface elevation, $\rho$ is the density, g is the gravitational acceleration, U, V and W are depth-averaged velocities in the x, y and z directions, respectively, q is the non-hydrostatic pressure at the bed surface, and S_x and S_y are source terms including the bottom frictions and viscous stresses, respectively [44,45,46]. The pressure is decomposed into hydrostatic (p_h) and non-hydrostatic (q) components as in [11,12,13,14,15,16,17,24,25,26,27,28,29,30,31,32,44,45,46] and Keller-box scheme [47] for vertical gradient approximation of the non-hydrostatic pressure in the formulation.

$$p=\rho g\left(\eta -z\right)+q=p_h+q$$

(5)

Equations (1)–(4) is solved with linear distributions assumed in the vertical direction for both the non-hydrostatic pressure and the vertical velocity.

2.2. Numerical Method and Least-Squares Finite-Element Method

A semi-implicit method similar to that in [11,12,13,14,15,16,17,24,25,26,27,28,29,30,31,32,44,45,46] is employed in this study. The depth-integrated non-hydrostatic SWM is solved semi-implicitly: (1) the provisional flow velocity is first implicitly solved using the hydrostatic shallow-water equations; (2) the non-hydrostatic pressure, which is implicitly obtained by ensuring a divergence-free velocity field, is used to correct the provisional velocity, and (3) finally the depth-integrated continuity equation is explicitly solved for the free-surface elevation to satisfy global mass conservation [11,12,13,14,15,16,17,24,25,26,27,28,29,30,31,32,44,45,46]. θ-method and least-squares finite-element method are used for time-stepping and spatial discretization, respectively.

After linearizing the nonlinear terms and θ-method time-stepping, Equations (6)–(8) become

$$\begin{array}{l}{\eta }^{n+1}+\theta \Delta t{\left[\left(\frac{\partial \tilde{U}}{\partial x}+\frac{\partial \tilde{V}}{\partial y}\right)\eta +\left(\tilde{U}\frac{\partial \eta }{\partial x}+\tilde{V}\frac{\partial \eta }{\partial y}\right)+\left(\frac{\partial \tilde{H}}{\partial x}U+\tilde{H}\frac{\partial U}{\partial x}\right)+\left(\frac{\partial \tilde{H}}{\partial y}V+\tilde{H}\frac{\partial V}{\partial y}\right)\right]}^{n+1}\\ ={\eta }^n+\theta \Delta t{\left[\left(\frac{\partial \tilde{U}}{\partial x}+\frac{\partial \tilde{V}}{\partial y}\right)h+\left(\tilde{U}\frac{\partial h}{\partial x}+\tilde{V}\frac{\partial h}{\partial y}\right)+\left(\frac{\partial \tilde{H}\tilde{U}}{\partial x}+\frac{\partial \tilde{H}\tilde{V}}{\partial y}\right)\right]}^n-\left(1-\theta \right)\Delta t{\left(\frac{\partial HU}{\partial x}+\frac{\partial HV}{\partial y}\right)}^n\end{array}$$

(6)

$$\begin{array}{l}{U}^{n+1}+\theta \Delta t{\left[\left(\tilde{U}\frac{\partial U}{\partial x}+\frac{\partial \tilde{U}}{\partial x}U+\tilde{V}\frac{\partial U}{\partial y}+\frac{\partial \tilde{U}}{\partial y}V\right)+g\frac{\partial \eta }{\partial x}\right]}^{n+1}+\theta \Delta t{S}_x^{n+1}\\ =U^n+\theta \Delta t{\left[\left(\tilde{U}\frac{\partial \tilde{U}}{\partial x}+\tilde{V}\frac{\partial \tilde{U}}{\partial y}\right)\right]}^n-\left(1-\theta \right)\Delta t{\left(U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}+g\frac{\partial \eta }{\partial x}\right)}^n+\left(1-\theta \right)\Delta t{S}_x^n\end{array}$$

(7)

$$\begin{array}{l}{V}^{n+1}+\theta \Delta t{\left[\left(\tilde{U}\frac{\partial V}{\partial x}+\frac{\partial \tilde{V}}{\partial x}U+\tilde{V}\frac{\partial V}{\partial y}+\frac{\partial \tilde{V}}{\partial y}V\right)+g\frac{\partial \eta }{\partial y}\right]}^{n+1}+\theta \Delta t{S}_y^{n+1}\\ =V^n+\theta \Delta t{\left[\left(\tilde{U}\frac{\partial \tilde{V}}{\partial x}+\tilde{V}\frac{\partial \tilde{V}}{\partial y}\right)\right]}^n-\left(1-\theta \right)\Delta t{\left(U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}+g\frac{\partial \eta }{\partial y}\right)}^n+\left(1-\theta \right)\Delta t{S}_y^n\end{array}$$

(8)

where superscript “n” and “n + 1” denote value of the current time step and new time step, respectively, with $t^{n+1}=t^n+\mathsf{\Delta }t$ and the time step size $\Delta t$; label “~” is annotated as the value of previous iteration or the value of “n” time step. By the interpolations of finite-element method, the unknowns $\underset{˜}{u}={\left\{\eta ,U,V\right\}}^T$ can be approximated as

$$\left\{\begin{array}{c}\eta \left(x,y\right)\\ U\left(x,y\right)\\ V\left(x,y\right)\end{array}\right\}=\left[\begin{array}{c}N\left(x,y\right)\\ N\left(x,y\right)\\ N\left(x,y\right)\end{array}\right]\left\{\begin{array}{c}\eta \\ U\\ V\end{array}\right\},$$

(9)

where N (x, y) is the spatial interpolation function. Substituting the approximations in Equation (9) into Equations (6)–(8), the residuals can be written as

$$\begin{array}{l}\left\{\begin{array}{c}{R}_{\eta }\\ R_U\\ R_V\end{array}\right\}=\left[\begin{array}{ccc}\left(1+\theta \Delta t\left(\frac{\partial \tilde{U}}{\partial x}+\frac{\partial \tilde{V}}{\partial y}\right)\right)N+\theta \Delta t\left(\tilde{U}\frac{\partial N}{\partial x}+\tilde{V}\frac{\partial N}{\partial y}\right)& \theta \Delta t\left(\frac{\partial \tilde{H}}{\partial x}N+\tilde{H}\frac{\partial N}{\partial x}\right)& \theta \Delta t\left(\frac{\partial \tilde{H}}{\partial y}N+\tilde{H}\frac{\partial N}{\partial y}\right)\\ \theta \Delta tg\frac{\partial N}{\partial x}& \left(1+\theta \Delta t\frac{\partial \tilde{U}}{\partial x}\right)N+\theta \Delta t\left(\tilde{U}\frac{\partial N}{\partial x}+\tilde{V}\frac{\partial N}{\partial y}\right)& \theta \Delta t\frac{\partial \tilde{U}}{\partial y}N\\ \theta \Delta tg\frac{\partial N}{\partial y}& \theta \Delta t\frac{\partial \tilde{V}}{\partial x}N& \left(1+\theta \Delta t\frac{\partial \tilde{V}}{\partial y}\right)N+\theta \Delta t\left(\tilde{U}\frac{\partial N}{\partial x}+\tilde{V}\frac{\partial N}{\partial y}\right)\end{array}\right]{\left\{\begin{array}{c}\eta \\ U\\ V\end{array}\right\}}^{n+1}\\ -{\left\{\begin{array}{c}\eta +\theta \Delta t\left[\left(\frac{\partial \tilde{U}}{\partial x}+\frac{\partial \tilde{V}}{\partial y}\right)h+\left(\tilde{U}\frac{\partial h}{\partial x}+\tilde{V}\frac{\partial h}{\partial y}\right)+\left(\frac{\partial \tilde{H}\tilde{U}}{\partial x}+\frac{\partial \tilde{H}\tilde{V}}{\partial y}\right)\right]-\left(1-\theta \right)\Delta t\left(\frac{\partial HU}{\partial x}+\frac{\partial HV}{\partial y}\right)\\ U+\theta \Delta t\left(\tilde{U}\frac{\partial \tilde{U}}{\partial x}+\tilde{V}\frac{\partial \tilde{U}}{\partial y}\right)-\left(1-\theta \right)\Delta t\left(U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}+g\frac{\partial \eta }{\partial x}\right)\\ V+\theta \Delta t\left(\tilde{U}\frac{\partial \tilde{V}}{\partial x}+\tilde{V}\frac{\partial \tilde{V}}{\partial y}\right)-\left(1-\theta \right)\Delta t\left(U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}+g\frac{\partial \eta }{\partial y}\right)\end{array}\right\}}^n\end{array}$$

(10)

According to the definition of the least squares method, by minimizing the residuals, it is expressed as

$$\min {\displaystyle {\int }_{\Omega }{\underset{˜}{R}}^2}d\Omega$$

(11)

where $\underset{˜}{R}={\left\{R_{\eta },R_U,R_V\right\}}^T$ are the residuals and Ω is the spatial domain. Equation (6) is equivalent to

$${\displaystyle {\int }_{\Omega }\frac{\partial \underset{˜}{R}}{\partial \underset{˜}{u}}\underset{˜}{R}d\Omega }$$

(12)

Details of the least-squares finite-element method can be found in related literatures [41,42].

Once the provisional flow velocity is obtained, a Poisson equation for the non-hydrostatic pressure can be formulated by the continuity equation with the provisional flow velocity, which is solved to ensure a divergence-free velocity field [11,12,13,14,15,16,17,24,25,26,27,28,29,30,31,32,44,45,46]. The obtained non-hydrostatic pressure is used to correct the provisional velocity

$$U^{c,n+1}=U^{n+1}+\Delta U^{n+1}=U^{n+1}-\frac{\Delta t}{2\rho }\left[\frac{\partial q^{n+1}}{\partial x}+\frac{q^{n+1}}{H^n}\left(\frac{\partial \eta }{\partial x}+\frac{\partial h}{\partial x}\right)\right]$$

(13)

$$V^{c,n+1}=V^{n+1}+\Delta V^{n+1}=V^{n+1}-\frac{\Delta t}{2\rho }\left[\frac{\partial q^{n+1}}{\partial y}+\frac{q^{n+1}}{H^n}\left(\frac{\partial \eta }{\partial y}+\frac{\partial h}{\partial y}\right)\right]$$

(14)

The superscript “c” denotes the correction. Equations (13) and (14) represents the final results of the flow velocity of the new time step. Free-surface elevation is then computed by the mass conservation Equation (15) with the corrected velocities

$${\eta }^{n+1}={\eta }^n+\Delta {\eta }^{n+1}={\eta }^n-\Delta t\left[\frac{\partial \left(H^n{U}^{c,n+1}\right)}{\partial x}+\frac{\partial \left(H^n{V}^{c,n+1}\right)}{\partial y}\right]$$

(15)

Details of the semi-implicit procedure can be found in [11,12,13,14,15,16,17,24,25,26,27,28,29,30,31,32,44,45,46].

3. Results

The non-hydrostatic SWM is verified with a solitary wave propagating in a channel in Section 3.1 and a sinusoidal wave passing over a submerged trapezoidal bar in Section 3.2, and applied to simulate von Karman vortex street in Section 3.3 and the ocean circulation of Dongsha Island water in Section 3.4.

3.1. Solitary Wave Propagation in a Channel

The solitary wave propagation is a nonlinear wave with finite amplitude. If the fluid is inviscid, the wave should maintain the shape and velocity during propagation process over a horizontal bottom. Solitary wave in a one-dimensional channel with a constant depth is used for model verification. An approximation of a solitary wave for the shallow-water equations is

$$\eta \left(x,t\right)=a\sec h^2\left(\sqrt{\frac{3a}{4H^3}}\left(x-x_0-ct\right)\right)$$

(16)

$$u\left(x,t\right)=\frac{\eta c}{\eta +h}\mathrm{and}w\left(x,t\right)=-z\frac{\partial u}{\partial x}$$

(17)

where a is the amplitude, h is the constant water depth, and $c=\sqrt{g\left(a+h\right)}$ is the wave propagation speed.

Simulation of the propagation of a solitary wave over a constant water depth H = 1.0 m with a = 0.1 m in a channel of length L = 80 m is performed. The channel is divided into 2000 uniform quadrilateral elements, i.e., $\mathsf{\Delta }x=\Delta y=0.04$ m and $\mathsf{\Delta }t=0.005$ s are used in computations. The corresponding Courant–Friedrichs–Lewy number CFL ≅ 0.4. Simulation starts with the wave profile and speed of Equations (16) and (17) at $x_0=10.0$ m and t = 0 s as the initial conditions.

Figure 1 shows the evolution of the solitary wave at t = 6 and 12 s. The profile of free-surface (η) and velocity (U) of the traveling solitary wave obtained with the hydrostatic SWM becomes distorted, the leading front becomes steeper, and the trailing front becomes smoother during the propagation progress. However, the profile of free-surface and velocity of the traveling solitary wave computed with the non-hydrostatic SWM remains close to that of the approximation—the wave shape is well conserved during wave propagation. There is a small reduction of wave amplitude and magnitude of velocity at the beginning of simulation because the initial conditions approximated by the analytical solution, Equations (16) and (17), is not the exact solution of the non-hydrostatic shallow-water equations, Equations (1)–(4). Similar observations were reported by [44,45].

Figure 2 shows the non-hydrostatic pressure (q) and vertical velocity (W) at t = 6 and 12 s obtained by the non-hydrostatic model, and the hydrostatic SWM does not have this information. It can be seen there are minor oscillatory dispersive trailing waves. Overall, computed results of non-hydrostatic SWM better resolve the salient characteristics of the propagating solitary wave than that of the hydrostatic SWM.

3.2. Propagation of Progressive Sinusoidal Waves over a Submerged Trapezoidal Bar

A physical experiment of progressive sinusoidal waves passing over a submerged trapezoidal bar has been conducted [48], as illustrated in Figure 3. It has been used to verify several non-hydrostatic models [29,44,45,46,49]. The objective is to test the capability of the non-hydrostatic models to simulate the relative strong interactions between nonlinear waves with the uneven bottom. The bottom topography is one of the major controlling factors, leading to significant reflection and transformation in wave field [38,50,51,52]. Besides, increased nonlinearity would further complicate their harmonic generation [48,53,54].

Figure 3 illustrates the geometry and experimental setup. The still water depth h = 0.4 m. Progressive sinusoidal waves with an amplitude a = 0.01 m and period T = 2 s are specified at the left inflow boundary, and they exit at the channel at right of the outflow boundary. Hydrostatic pressure assumption for shallow-water flows is generally valid for water depth parameter kh < 0.31. The values of kh in region I and II are 0.68 and 0.32, respectively. The value of kh ranges from 6 to 10 in region III based on a spectral analysis [54].

The computational domain is discretized by a set of uniform of 875 quadrilaterals and 1752 nodes, equivalent to grid size $\Delta x$ = $\mathsf{\Delta }y$ = 0.04 m, and a time step size of $\Delta t$ = 0.004 s is chosen. Time history of computed free-surface elevation is plotted and compared with the experimental data of [48] and prediction of the five $\delta$-layers non-hydrostatic model of [29,49] at six different wave gauge locations (gauge 2–7) in Figure 4.

When the wave trains propagate onto the submerged step, the wave amplitude and steepness increase due to shoaling. Though model predicted shoaling coefficient (equivalent to the relative wave height $\max \left|\eta \right|/a={\left[2{\cos \mathrm{h}}^{2}kh/\left(\sin \mathrm{h}2kh+2kh\right)\right]}^{1/2}$) 1.23 at gauge 3 agrees well with 1.34 by the Airy linear wave analytic solution [55], the model underestimates steeper shoaling waves on the upward slope at gauge 2 and 3 compared with the measurement data. Apparently, wave transformation here is a nonlinear process, not linear. The model captures the main features of the riding waves over the bar, wave trains travel with the correct phase speed, and higher harmonics are generated and secondary waves are developed, but are not obvious, at gauge 4 and 5. However, the predicted wave peak value is much smaller than the measured value and computed value of [29,49]. Behind the step, the wave amplitude and steepness decrease at station 6 and 7. The nonlinearity of wave becomes weaker, yielding the release of the generated higher harmonic waves where wave dispersion is important. There are clearly visible discrepancies between prediction of non-hydrostatic SWM with the measurement data [48] and prediction of 5$\delta$-layers non-hydrostatic SWM [29,49]. These discrepancies mainly due to kh ranging from 6 to 10 in region III is out of the applicable range of the current one-layer depth-integrated non-hydrostatic SWM. Similar underestimates and discrepancies are reported in [44,45]. One-layer vertical integral non-hydrostatic SWM is not able to accurately resolve the details of vertical flow structures where nonlinearity and dispersion are important. Moreover, predictions with the hydrostatic SWM, neglecting the dispersion terms, is totally unrealistic and different from the predictions with measured data [48] and 5$\delta$-layers non-hydrostatic SWM [29,49].

3.3. von Karmann Vortex Street

The wake flow behind a bluff obstacle has been a subject of interest to engineers as well as to scientists for many years. The periodic alternate shedding of counter-rotating vortex pairs downstream the bluff obstacle is the so-called von Karmann vortex street. It is a classic fluid dynamic problem resulted from flow passing an obstacle at a certain range of Reynold’s number, i.e., 46 < Re < 1000 [56,57,58]. Bluff body wakes have significant different characteristic flow properties depending on the Reynolds number Re = $U_r{L}_r/{\nu }_r$, where $U_r$ is the characteristic velocity, $L_r$ is the characteristic length, and ${\nu }_r$ is the eddy viscosity, respectively.

Both hydrostatic and non-hydrostatic shallow-water model are used to simulate the von Karmann Vortex shedding due to passing of uniform flow over a square obstacle. The length and width of the study domain are 400 and 200 m, respectively. A square obstacle with length and width equal to 20 m is located in the center of the mainstream flow. The still water depth is constant h = 4 m. Flow enters the study domain from the left boundary with a constant speed $U_r$ = 0.2 m/s and exits the study domain from the right open boundary. Free slip boundary condition is employed for both top and bottom boundaries. The resulting Reynolds number Re = $U_r{L}_r/{\nu }_r$ = 400, where $U_r$ = 0.2 m/s is the characteristic velocity, $L_r$ = 20 m is the characteristic length (square obstacle dimension), and ${\nu }_r$ = 0.01 m²/s is the eddy viscosity, respectively.

Figure 5a depicts the study domain and computational meshes, consisting of 20,996 nodes and 20,754 quadrilaterals. Fine meshes are used near the square obstacle region, Figure 5b, where flow and wave field are expected to change more significantly than other areas. Simulation starts with the still water, i.e., fluid is motionless initially. It takes a short period of transition for the flow to develop and reach a periodic state.

Figure 6 shows the comparison of time history of velocity components of point P, located at (25, 15) and depicted in Figure 5b, by hydrostatic and non-hydrostatic SWM. As can be seen, the difference of predicted velocity components by hydrostatic and non-hydrostatic SWM is not significant—the amplitude of velocity components remains the same and phase shift increases in the long period of computations. The period (T) of oscillation of the von Karman vortex street is 480 s for hydrostatic SWM and 475 s for non-hydrostatic SWM, and the corresponding Strohal number St = L_r/U_rT ≈ 0.208 and 0.210 for hydrostatic and non-hydrostatic SWM, respectively. This value is close to 0.218 obtained by the formula $St=0.2698-1.0271/\sqrt{\mathrm{Re}}$ given by [59].

Since the difference of prediction of hydrostatic and non-hydrostatic SWM is not significant and flow fields of the two models are visually indistinguishable, only the simulation results of hydrostatic SWM are shown. Figure 7 and Figure 8 show the contours of velocity components and streamlines at t = T/2 and T. A well-organized von Karman vortex street is reproduced.

3.4. Ocean Circulation of Dongsha Water

The Dongsha Islands, also known as Pratas Islands, are three atolls in the north of the South China Sea. They consist of one island, two coral reefs and two banks as well as lagoon, and channels. Figure 9 depicts the location of Dongsha Atolls and the surrounding water. There are gaps in the northwest and southwest of the atolls. The Dongsha Island, located in the middle of gaps, forms the northern and southern channels which leads (exits) water into (from) the atolls at northern and southern sides of Dongsha Island, respectively.

There is a sand spit which is commonly referred to as the “swaying dragon tail”. The sand spit takes shape from the process of alongshore drift and sedimentation by nearshore currents at the headland. The phenomenon of swaying dragon tail is influenced by waves and tides induced by the northeast monsoons and southwest monsoons. The sediment gradually drifts from north to south in the winter and from south to north in the summer. By the in-situ observation, it has been shown that the swaying dragon tail suffers the inappropriate coastal protections such as the constructions of armor units, groins and detached breakwaters around the island, leading to the vanishing of some special landscapes.

Coastal erosion is a serious problem in Dongsha Island by natural factors (e.g., sea level rise) as well as human being actions. Recent measurements [62] of the Dongsha Atolls show that the water depth is about 2–22 m and the averaged velocity is 0.06 m/s and 0.24 m/s in the internal atoll and external atoll, respectively. During the flood tides, the counterclockwise current mainly exchanges water masses through the northern channel to the internal atoll. At the same time, the southern channel emerges a weaker current and presents an uneven spatial distribution of speed. An opposite-direction current exits near the southern channel and northern channel in internal atoll during flood tides. During the ebb tides, the mainstream current exits from northern channel, and the southern channel produces the vortexes repeatedly.

POM (Princeton Ocean Model) is widely used by oceanic and atmospheric models. In the early 1980s, the model was used primarily to coarse-resolution simulations of the large-scale ocean circulation. With the advancement of computer hardwares and numerical algorithms, the model is able to handle high-resolution coastal ocean processes [63]. POM has been applied to simulate the ocean circulation and transportation of substances of Dongsha waters [64]. However, the flow field of the internal atoll is not considered in the model due to the coarse mesh resolution used. A high resolution two-dimensional non-hydrostatic shallow-water model is applied to study the hydrodynamics of the Dongsha water. The computational domain length and width are 152 km and 92 km, respectively. Figure 10 plots the computational meshes, total 35,923 nodes and 17,788 elements. The complex topography and coastline is approximated using the smaller meshes to achieve the accuracy, as shown in Figure 10a. The high resolution of meshes is applied around the atoll region where large change of flow field is expected. In contrast, the relative large meshes are used for smooth changing flow field areas. The interpolation method of IDW (inverse distance weighted) is used to interpolate water depth from ETOPO1 [60] to non-hydrostatic SWM, shown in Figure 10b.

The harmonic analysis is similar to the spectral analysis of investigation by [62]. The tide is mainly composed of diurnal tides K₁ and O₁ and semi-diurnal tides M₂. The period of diurnal tide K₁ is used in the modeling. Boundary conditions used are no flux across the north and south boundaries, open boundary condition on the east boundary, and a K₁ tide is specified on the west side in the Dongsha ocean circulation modeling. Simulation of flow starts from the static state, and it takes a period of time to reach a periodic state. Parameters used for simulations of Dongsha ocean circulations are listed in

Table 1. Parameters used in the Dongsha ocean circulations modeling.

Length of Computed Domain (L).	Width of Computed Domain (W)	Period of K1 Tide (T)	Mesh Size (Δx)	Time Step (Δt)
152 km	92 km	86,160 s	120–5000 m	100 s

.

Figure 11a,b shows the measured flow field at ebb tides and flood tides [62]. Computed flow fields at ebb and flood tides by non-hydrostatic SWM are similar to the flow field of measurements, as shown in Figure 11c,d. The concentrated flow mainly exits from northern channel and a clockwise moving vortex is observed along the northern edge of the channel at ebb tides. The concentrated flow mainly enters from the northern channel and a counterclockwise moving vortex is observed along the edge of the northern channel at flood tides. The flow field presents an opposite pattern in the southern channel, and the flow is weaker in the southern channel compared with that in the northern channel due to the influence of topography in the southern channel at flood tides.

A computational mesh of POM, total of 8447 nodes, 14,999 triangular elements, and 20 non-uniform vertical layers, has been applied to Dongsha water. Bathymetry of the atoll is assumed to be shallow for computational stability. Figure 12a–c illustrates the distribution of flow field at ebb tides and flood tides by POM, respectively. The flow moves from west to east in ebb tides. The velocities are large at northern and southern channels and around the Dongsha Island. In contrast, in the flood tides, the flow moves from east to west. Figure 12a–d shows the comparison of predicted flow field of POM and non-hydrostatic SWM in ebb tides and flood tides, respectively.

Periodically moving vortexes persist around the upper and lower edge of the northern and southern channels of Dongsha Atoll. The mainstream flow exits the interior atoll from the northern channel at ebb tides. In the meantime, there is a pair of moving vortices—a clockwise vortex in the north edge of the northern channel and a counterclockwise vortex in the south edge of the southern channel. The mainstream flows enter the interior of atoll from the northern channel at flood tides, and the moving vortices change the direction of rotations comparing that of ebb tides. The clockwise moving vortex in the southern channel becomes weaker because of the effect of terrain at flood tides. General features of simulated flow field compared reasonably similar to that of the flow field of measurements. Value of computed speed 0.1–0.3 m/s by non-hydrostatic SWM is smaller than 0.8–1.0 m/s by the POM. The discrepancy is mainly attributed to the modification of atoll and the top-layer flow field of POM compared with the depth-averaged flow field of non-hydrostatic SWM.

4. Discussion

Hydrostatic pressure assumption for shallow-water flows is generally valid for water depth parameter kh < 0.31. Hydrostatic SWM is efficient for modeling shallow water flows where variation of vertical scale is much smaller than the variation of horizontal scale. To relax the constraints of hydrostatic pressure assumption and non-dispersion, a non-hydrostatic shallow-water model is developed for weakly dispersive waves.

Both wave and flow driven fluid flows, including solitary wave propagation in a channel, propagation of progressive sinusoidal waves over a submerged bar, von Karmann vortex street, and ocean circulations are numerically investigated. It is found that the difference between predictions of hydrostatic and non-hydrostatic SWM is insignificant for flow-induced fluid flows, where dispersion plays a minor role in these situations. However, for wave-induced fluid flows, difference between predictions of hydrostatic and non-hydrostatic SWM can become pronounced, especially for nonlinear waves traveling from deep water to shallow water and encountering the uneven bathymetry.

Hydrostatic SWM cannot accurately model the solitary wave propagation in a constant water depth channel and progressive sinusoidal waves propagation over a submerged trapezoidal bar, where both dispersion and nonlinearity play an important role. Non-hydrostatic SWM can accurately predict the solitary wave propagation, but cannot model well the progressive sinusoidal wave propagation over a submerged trapezoidal bar, because the water depth parameter kh ranging from 6 to 10 behind the bar is out of applicable range of the present one-layer depth-integrated non-hydrostatic SWM. Therefore, to accurately simulate the nonlinear dispersive free-surface flows, a multi-layer, non-hydrostatic shallow-water model following [26,27] is suggested for the future study.

Hydrostatic SWM (HSWM) is computationally efficient for modeling shallow water flows than the non-hydrostatic SWM (NHSWM) modeling. In general, NHSWM takes 1.5 to 3.0 more cpu time than HSWM does for the same simulation. The most cpu consuming part of NHSWM computations is to solve the Poisson-type equation for the non-hydrostatic pressure. It takes about 40–70% of the cpu time of the computations, depending on the complexity of wave/flow fields, number of degree-of-freedoms (nodes and elements), treatment of the boundary conditions, and convergence criterion, etc. Moreover, the increase in the percentage of cpu time for solving the non-hydrostatic pressure equation is larger than linear proportion when increasing the number of nodes and elements.

Since the least-squares formulation is employed, the resulting system of equations of Poisson-type non-hydrostatic pressure is symmetric and positive-definite. Therefore, Jacobi precondition with element-by-element conjugate gradient method [43] is used to solve the resulting system of equations. Furthermore, numerical solution of the non-hydrostatic pressure is prone to numerical oscillations. To suppress the oscillations, either a finer mesh resolution or smoothing is employed. In this study, the Laplacian smoothing is used.

For applications involving short period waves, abruptly changing bed topographies, and stratification due to strong density gradients, the hydrostatic assumption is invalid. To accurately simulate the nonlinear dispersive free-surface flows, a multi-layer non-hydrostatic shallow-water model following [26,27] is suggested for the future study.

5. Conclusions

A non-hydrostatic shallow-water model is developed and applied first to three benchmark tests and then to the ocean circulation surrounding Dongsha Atolls of the South China Sea. Least-Squares Finite-Element Method (LSFEM) is used to solve the non-hydrostatic shallow-water equations. There are many advantages of LSFEM compared with popular Galerkin and other methods, such as a single approximating space can be used for all variables and the choice of approximating spaces is not subject to the Ladyzhenskaya–Babuska–Brezzi (LBB) condition. Therefore, mixed interpolation formulation, collocation grids and convection scheme are not necessary. More importantly, the resulting system of equations of LSFEM is symmetric, positive-definite (SPD), and therefore, can be used with an efficient element-to-element conjugate gradient method.

Difference between predictions of hydrostatic and non-hydrostatic SWM for von Karmann vortex street and ocean circulation of Dongsha water is found insignificant. However, computations with the non-hydrostatic SWM require much more computer resources and cpu time compared with computations with the hydrostatic SWM. For instance, it takes about twice cpu time to model the von Karmann vortex street with non-hydrostatic SWM than with hydrostatic SWM, and about 2.5 times cpu time to model the ocean circulation of Dongsha water with non-hydrostatic SWM than with hydrostatic SWM.