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

: A two-dimensional non-hydrostatic shallow-water model for weakly dispersive waves is developed using the least-squares ﬁnite-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 ﬁnite-element method. In the predictor step, hydrostatic pressure at the previous step is used as an initial guess and an intermediate velocity ﬁeld is calculated. In the corrector step, a Poisson equation for the non-hydrostatic pressure is solved and the ﬁnal velocity and free-surface elevation is corrected for the new time step. The non-hydrostatic shallow-water model is veriﬁed and applied to both wave and ﬂow driven ﬂuid ﬂows, 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 e ﬃ cient and accurate for shallow water ﬂows. 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 a ﬀ ect the velocity ﬁeld and free-surface signiﬁcantly in case where nonlinearity and dispersion are important during the course of wave propagation.


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

Shallow-Water Equations
The two-dimensional depth-averaged non-hydrostatic shallow-water equations are derived from the depth-averaged Navier-Stokes equations and read.
Equations (1)-(4) is solved with linear distributions assumed in the vertical direction for both the non-hydrostatic pressure and the vertical velocity.
After linearizing the nonlinear terms and θ-method time-stepping, Equations (6)- (8) become where superscript "n" and "n + 1" denote value of the current time step and new time step, respectively, with t n+1 = t n + ∆t and the time step size ∆t; label "~" is annotated as the value of previous iteration or Water 2020, 12, 2832 4 of 17 the value of "n" time step. By the interpolations of finite-element method, the unknowns u = η, U, V T can be approximated as 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 According to the definition of the least squares method, by minimizing the residuals, it is expressed as are the residuals and Ω is the spatial domain. Equation (6) is equivalent to 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 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 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].

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.

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 where a is the amplitude, h is the constant water depth, and c = g(a + h) 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., ∆x = ∆y = 0.04 m and ∆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].
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.

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 where a is the amplitude, h is the constant water depth, and    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, Equation (16) and (17), is not the exact solution of the non-hydrostatic shallowwater 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. 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 nonhydrostatic SWM better resolve the salient characteristics of the propagating solitary wave than that of the hydrostatic SWM.

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 nonhydrostatic 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].

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].

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 nonhydrostatic 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].   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 ∆x = ∆y = 0.04 m, and a time step size of ∆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 δ-layers non-hydrostatic model of [29,49] at six different wave gauge locations (gauge 2-7) in Figure 4. The computational domain is discretized by a set of uniform of 875 quadrilaterals and 1752 nodes, equivalent to grid size x  =Δ = 0.04 m, and a time step size of 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 -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 heightmax| |/ = [2cosh ℎ/(sinh2 ℎ + 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]. 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 η /a = 2cos h 2 kh/(sin h2kh + 2kh) 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δ-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δ-layers non-hydrostatic SWM [29,49].

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 /ν r , where U r is the characteristic velocity, L r is the characteristic length, and ν r is the eddy viscosity, respectively.
Both 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 -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 - layers nonhydrostatic SWM [29,49].

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 = / , where is the characteristic velocity, is the characteristic length, and r  is the eddy viscosity, respectively.     Re given by [59]. 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 = Lr/UrT  0.208 and 0.210 for hydrostatic and non-hydrostatic SWM, respectively. This value is close to 0.218 obtained by the formula = 0.2698 − 1.0271/√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. Figures 7 and 8 show the contours of velocity components and streamlines at t = T/2 and T. A well-organized von Karman vortex street is reproduced. 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 = Lr/UrT  0.208 and 0.210 for hydrostatic and non-hydrostatic SWM, respectively. This value is close to 0.218 obtained by the formula = 0.2698 − 1.0271/√Re given by [59].

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

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.

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 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. 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 oppositedirection 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 K1 and O1 and semi-diurnal tides M2. The period of diurnal tide K1 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 K1 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. The harmonic analysis is similar to the spectral analysis of investigation by [62]. The tide is mainly composed of diurnal tides K 1 and O 1 and semi-diurnal tides M 2 . The period of diurnal tide K 1 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 1 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. Table 1. Parameters used in the Dongsha ocean circulations modeling.

Length of Computed Domain (L). Width of Computed Domain (W) Period of K 1 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.
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 nonuniform 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 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.
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 toplayer flow field of POM compared with the depth-averaged flow field of non-hydrostatic SWM.

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.

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.