the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# A nonhydrostatic oceanic regional model, ORCTM v1, for internal solitary wave simulation

### Hao Huang

### Pengyang Song

### Shi Qiu

### Jiaqi Guo

### Xueen Chen

The Oceanic Regional Circulation and Tide Model (ORCTM), including a nonhydrostatic dynamics module which can numerically reproduce internal solitary wave (ISW) dynamics, is presented in this paper. The performance of a baroclinic tidal simulation is also examined in regional modeling with open boundary conditions.

The model control equations are characterized by three-dimensional and
fully nonlinear forms considering incompressible Boussinesq fluid in
*Z* coordinates. The pressure field is decomposed into the surface,
hydrostatic, and nonhydrostatic components on the orthogonal curvilinear
Arakawa-C grid. The nonhydrostatic pressure determined by the intermediate
velocity divergence field is obtained via solving a three-dimensional
Poisson equation based on a pressure correction method. Model validation
experiments for ISW simulations with the topographic change in the
two-layer and continuously stratified ocean demonstrate that ORCTM has a
considerable capacity for reproducing the life cycle of internal solitary
wave evolution and tide–topography interactions.

Internal wave (also called internal gravity waves) activities have been observed frequently across the stratified ocean and play a significant role in the multiscale energy cascade (Mtfller, 1976). Observations reveal that internal waves, especially high-frequency internal solitary waves, could contain significant potential energy with strong vertical shear, mixing, and wave breaking, leading to a dramatic change in the currents and density structures (Ramp et al., 2004; Vlasenko et al., 2010; Huang et al., 2016), violent overturning bringing sediment and nutrients from the seafloor to the surface (Wang et al., 2007), and even damage to some underwater vehicles (Duda et al., 2006) and deep-water drilling projects (Osborne et al., 1978). Basically, astronomical tides passing abrupt topography can cause the generation of baroclinic tides (also called internal tides, hereafter ITs) with multi-modal structures then capable of propagation, disintegration, and dissipation in the ocean (Vlasenko et al., 2005, 2010). The low mode of baroclinic tides can travel thousands of kilometers with long horizontal wavelengths of tens of kilometers (Baines, 1982). Furthermore, the inclusion of nonlinear and nonhydrostatic effects permits the evolution of nonlinear internal waves (hereafter NIWs) and even internal solitary waves (hereafter ISWs) derived from the steepening of low-mode internal tides as the consequence of the ever-changing terrain and background stratification (Gerkema and Zimmerman, 1995; Legg and Adroft, 2003).

Numerical ocean models are one of the most effective tools to study internal waves compared to theoretical methods, in situ observations, and laboratory investigations. Ocean models with the hydrostatic balance approximation have been used to explore the regional circulation and tide dynamics across the stratified ocean. The hydrostatic balance manages to take the large scales to mesoscales into consideration due to the fairly high accuracy (Marshall et al., 1997b; Chen et al., 2003; Shchepetkin and McWilliams, 2005; Ko et al., 2008). However, In the hydrostatic balance scheme, omitting some essential terms in the vertical momentum equation results in the inapplicability of nonhydrostatic dynamics (Lai et al., 2010). For example, the subsequent steepening of the internal tides and the later high-frequency nonlinear ISW formation cannot be depicted by a hydrostatic modeling whereby only internal jumps are formed but no soliton appears (Li, 2014) because the strong vertical current with its order of magnitude equals the horizontal one via the scale analysis method (Marshall et al., 1997a). In other words, the three-dimensional Navier–Stokers equations should be considered thoroughly. It is indispensable for simulating the nonlinear and large-amplitude ISWs to develop a nonhydrostatic ocean model in consideration of nonhydrostatic dynamics.

A robust ocean model with nonhydrostatic dynamics realizations should satisfy at least two requirements synchronously: (1) high enough accuracy of mesoscale to large-scale simulations must be guaranteed, such as large-scale wind-induced circulation and mesoscale eddies reconstructed and mainly influenced by the hydrostatic balance. (2) Meanwhile, they should be concerned with small scales to mesoscales with higher spatial and temporal resolution that are resolved finely under the nonhydrostatic balance; for instance, there is a simulation able to describe the cradle-to-grave process for tide–topography interactions, the dispersive effects and nonlinear steepening of baroclinic tides, and the breaking and dissipation of strong nonlinear ISWs. The nonhydrostatic simulation can apply to the small to large scales across the stratified ocean simultaneously, which is identified as one of the main directions for research and development of the nonhydrostatic ocean model. In reality, there have been some nonhydrostatic ocean models or ones considering nonhydrostatic dynamics coming out in the past decades, such as MITgcm (Marshall et al., 1997a, b, 1998), SUNTANS (Fringer et al., 2006), and ROMS (Kanarska et al., 2007). All of the above have been used to realize a series of two- or three-dimensional nonhydrostatic numerical studies, including the instability of small-scale flows in laboratory experiments (Lai et al., 2010; Li et al., 2022), internal solitary waves in the continental shelves (Vlasenko et al., 2010; Zeng et al., 2019), and the hydraulic lee wave around the seamount (Kanarska et al., 2007; Liu et al., 2016). Nevertheless, the primary reason why there is still no widespread use for the nonhydrostatic ocean model is that the nonhydrostatic solution to an extensive sparse linear equation is too demanding to solve directly for the 3-D oceanic environment. That usually demands large iteration times, fast convergence speed, and large PC storage. For this reason, Ai and Ding (2016) employed a novel model grid arrangement to render the sparse linear equation discretized form simpler to solve where the bottom-fitted coordinate ensures the homogeneous boundary condition. Moreover, numerical errors can be avoided via the immersed boundary method to treat uneven bottoms in the calculation of the baroclinic pressure force (Ai et al., 2021). Generally, whether the boundary conditions are matched with the whole nonhydrostatic algorithm can shape the performance of complex nonhydrostatic dynamics in the regional ocean model. In addition, the different kinds of sub-grid parameterization schemes have a profound impact on the model performance with a necessity for appropriate ones to be assessed, and most of these model codes are seldom shared or open-source. Supposing we develop a nonhydrostatic ocean model based on an original hydrostatic framework model. In that case, the nonhydrostatic dynamics module should involve a complete vertical momentum equation. Some terms associated with the vertical velocity are required to be complemented simultaneously in the other equation. Also, based on the idea of the fractional step method (Press et al., 1988; Armfield and Street, 2002), the total pressure is to be decomposed into hydrostatic and nonhydrostatic components (Marshall et al., 1997a; Lai et al., 2010). The former corresponds to the result of hydrostatic balance, and the divergence for intermediate velocity limits the latter to correct the local velocity fields, called the “pressure correction” method (Stansby and Zhou, 1998; Fringer et al., 2006; Kanarska et al., 2007; Lai et al., 2010). With these methods, the nonhydrostatic dynamics simulation can be economically fulfilled comparatively in harmony with the original physical framework as an extension of the hydrostatic ocean model.

In this context, we have implemented the nonhydrostatic dynamic algorithm into the Oceanic Regional Circulation and Tide Model (hereafter ORCTM) and demonstrated its capability and performance of reproducing the life cycle of nonlinear internal solitary waves in distinct hydrodynamic environments. The rest of the paper is organized as follows. In Sect. 2, the basic framework of ORCTM including control equations, open boundary conditions, and nonhydrostatic algorithms is described. In Sect. 3, a series of numerical validation experiment results is presented, aiming at the simulation of the overall processes of internal solitary waves. In the last section, we have some further discussions and come to conclusions.

The Max Planck Institute Ocean Model (MPI-OM) is a global ocean circulation
and tide model based on the ocean primitive equations discretized on the
orthogonal curvilinear Arakawa-C grid with hydrostatic balance approximation
(Marsland et al., 2003; Chen et al., 2005). Rooted from MPI-OM, in this
paper, an oceanic regional circulation and tide model (ORCTM) has been
developed to realize the simulation for nonhydrostatic internal solitary
wave modeling, which will be referred to hereafter as ORCTM version 1.0.
The *z*-level grid applied has the partial filled cell capability to adjust
the distance of the vertical grid on the seabed for fitting into the realistic
terrain, and the tidal forcing flow can be implemented via a relaxation
scheme at the open boundary with an area of sponge layers. It is under the
laws of the Boussinesq, rotating, and fully nonlinear Navier–Stokes fluid
that ORCTM can be used to reproduce and explore the nonhydrostatic dynamics
such as large-amplitude ISWs, nonlinear tidal internal waves, and
downwelling and upwelling processes of real oceans.

## 2.1 Control equations

The three-dimensional ocean primitive control equations involve the momentum, continuity, potential temperature, salinity, and density equations given as follows.

In the local Cartesian framework of reference on the rotating Earth for a
geophysical flow, *t* is the time; $\partial /\partial t$ is
the time partial derivative; *x*,*y*, and *z* axes are directed eastward, northward,
and upward, respectively; the horizontal velocity vector is
${\mathit{u}}_{h}=\left(u,v\right)$; and *w* is the
vertical velocity. With the linearized kinematic boundary condition and the
freshwater forcing term *Q*_{ς} from the evaporation and
precipitation (Marsland et al., 2003), the free surface elevation equation
can be proposed as follows.

*ς*
is the change in the free surface elevation; *P*, *θ*, and
*S* are pressure, potential temperature, and salinity; and *ρ*_{c} is the
reference density of seawater. The first and second Coriolis parameters
are *f*=2Ωsin *φ* and
$\stackrel{\mathrm{\u0303}}{f}=\mathrm{2}\mathrm{\Omega}\mathrm{cos}\mathit{\phi}$, where Ω is the
rotational angular speed and *φ* is the geographic latitude.
∇_{H} is the horizontal divergence operator; *Q*_{θ}
and *Q*_{s} are source or sink terms about potential temperature and
salinity. The equation of seawater state is the polynomial form for the
density *ρ* advocated by the Joint Panel on Oceanographic Tables and
Standards (Fofonoff and Millard, 1983). The additional forcing term vector
$\mathit{F}=\left({F}_{x},{F}_{y},{F}_{z}\right)$ can consider tidal
potential forcing. The horizontal eddy viscosity vector is ${\mathit{F}}_{\mathrm{H}}=\left({F}_{\mathrm{H}x},{F}_{\mathrm{H}y},{F}_{\mathrm{H}z}\right)$
described with the scale-dependent biharmonic formulation (Wolff et al.,
1997; Marsland et al., 2003), and the horizontal diffusivity terms of
temperature and salinity are *F*_{Hθ} and *F*_{HS}, supporting the
harmonic forms. In addition, the vertical eddy viscosity vector is ${\mathbf{F}}_{\mathrm{V}}=\left({F}_{\mathrm{V}x},{F}_{\mathrm{V}y},{F}_{\mathrm{Vz}}\right)$ and
eddy diffusivity terms are *F*_{Vθ} and *F*_{VS}. Here, the vertical
eddy turbulent frictions are specified to depend on the Richardson number
*R**i* via the modified PP81 parameterization scheme (Pacanowski and
Philander, 1981). The viscous terms above are expressed as follows.

Here, $\mathrm{\Delta}={\mathrm{\nabla}}_{\mathrm{h}}\cdot {\mathrm{\nabla}}_{\mathrm{h}}$ is the
horizontal Laplace operator, *B*_{H} and *D*_{H} are parameterized with the
horizontal grid resolution, and *N*(z) is the buoyancy frequency.
${A}_{\mathrm{V}}^{n+\mathrm{1}}$ and ${D}_{\mathrm{V}}^{n+\mathrm{1}}$ are updated on Eqs. (11) and (12) with
the time relaxation coefficient *λ* at every time step. Apart from
the background viscous coefficients *A*_{b} and *D*_{b} due to internal
wave breaking, the modified PP81 scheme also considers the wind-induced
turbulent coefficients *A*_{w} and *D*_{w} associated with the local mixed
layer depth and 10 m wind speed (Marsland et al., 2003). Here, the constant
number *α* is set to be 5. And the adjustable parameters *A*_{V0} and
*D*_{V0} can be determined by estimating energy flux at every grid point. As
for the boundary condition, the slip conditions are specified at surface and
bottom boundaries where the wind stress *τ*_{w} is based on the model input, and the bottom drags
*τ*_{b} are described by
linear and quadratic functions (Arbic and Scott, 2008). The top and bottom
boundary conditions can be written as

where *γ* and *C*_{d} are the bottom friction and drag coefficients
representing the linear and quadratic relation expressions, respectively.

## 2.2 Settings of open boundary conditions

It is fundamental for the regional model to be configured by an open boundary condition that avoids reflection waves effectively so that the outward waves can freely flow through the boundaries. Meanwhile, external inputs such as tidal waves can stably force the model domain through the boundaries, satisfying the need for consistency in hydrodynamics and computational mathematics. Here, we use the relaxation boundary conditions with sponge layers by consulting Zhang et al. (2011) that can dampen the reflection waves back into the interior domain and avoid the sharp gradients of water properties caused by the prescribed values on the boundaries. Specifically, we add a relaxation term $M\left(x,y,z,t\right)$ formularized with the exponential function in the specified sponge zones. At each time step, the model variables are updated with an explicit scheme expressed as follows.

In Eqs. (15) and (16), *m*_{b} is the boundary value of requisite
model variables including velocity, potential temperature, and salinity; *m* is the corresponding relaxation result in the interiors; *m*^{∗} is
the intermediate variable; *r* is the distance from the boundary; and Δ*t* is the model time step. Here, it should be noted that *τ* and *L*_{sp} are artificially prescribed adjustment parameters referring to
the timescale coefficient and the thickness of the sponge relaxation
layers. The model target variables over the sponge layer will relax
exponentially to the boundary values through the relaxation term, with
relaxation modulated by *τ* and *L*_{sp} in the exponential
shape. To restrain the reflection of outflow current, *τ* and *L*_{sp} need to be determined in advance via estimating the energy flux
of internal signals through the boundaries. This open boundary relaxation
condition is suitable for the numerical study of the large-amplitude ISWs so
that the outward strong, nonlinear, and nonhydrostatic wave and current
signals will dampen gradually.

## 2.3 Implementation of nonhydrostatic algorithms

According to the momentum equations (Eqs. 1–3), the total pressure *P* consists
of sea surface pressure *p*_{s}, hydrostatic pressure *p*_{h}, and
nonhydrostatic pressure *p*_{nh} given as follows.

It is negligible for the change in the sea surface pressure term *p*_{s} to impact the water column if the external atmospheric forcing is
excluded. Hydrostatic pressure *p*_{h} can be calculated from the
hydrostatic balance equation (Eq. 18), and the vertical momentum equation (Eq. 3) at
this stage becomes

where the left term refers to the local time rate of change, and the right term is the sum of the other forces without the additional forcing term vector. Compared with Eq. (18), the vertical momentum equation (Eq. 19) can also be called the nonhydrostatic balance equation. Furthermore, with the idea of the fractional step method (Press et al., 1988; Kanarska et al., 2007), the intermediate velocity field $\stackrel{\mathrm{\u0303}}{\mathit{u}}=\left(\stackrel{\mathrm{\u0303}}{u},\stackrel{\mathrm{\u0303}}{v},\stackrel{\mathrm{\u0303}}{w}\right)$ will be updated via the nonhydrostatic pressure ${p}_{\mathrm{nh}}^{n}$ gradients, which can be obtained via Eqs. (20) to (22) discretized as follows.

Here, the superscript *n* means the current time step and the vector $\mathit{G}=\left({G}_{x},{G}_{y},{G}_{z}\right)$ represents
the sum of the advection term, Coriolis term, eddy viscosity term, and
hydrostatic pressure gradient term. The discretized partial equations (Eqs. 23–25) are subsequently established under the relationship between the
nonhydrostatic pressure perturbation ${{p}^{\prime}}_{\mathrm{nh}}$ gradients and the next
time step *n*+1 velocity field. Then nonhydrostatic pressure at the next time
step is defined as Eq. (26) in light of the pressure
correction method. To acquire nonhydrostatic pressure perturbation the
continuity equation (Eq. 4) needs to be substituted into Eqs. (23) to (25) to
eliminate the following time step *n*+1 velocity field with the
three-dimensional Poisson equation (Eq. 27), which demonstrates that the
nonhydrostatic pressure depends on the vanishing of the divergence-free
velocity fields.

The Poisson equation (Eq. 27) can be discretized into a linear matrix equation (Eq. 28); the right-hand-side **B** is determined by the divergence of
the intermediate velocity field. The adjoint matrix **A**
represents the discrete three-dimensional Laplacian operator with a size of
the number of model cells. Their specific discrete processes are introduced
in Appendix A.

The proper boundary conditions need to be given to solve this Poisson
equation (Eq. 27). Here, the homogeneous Neumann boundary condition at the solid
boundaries, also called the zero-gradient condition (Eq. 29), is used with good
compatibility with the no-flux normal to slope, where ** n** is the
normal unit vector (Marshall et al., 1997a). We assume that nonhydrostatic
dynamic processes are weak enough at the sea surface and open boundaries. In
other words, the input signals through the boundaries are dominantly
hydrostatic with nonhydrostatic pressure perturbation close to zero. The
nonhydrostatic dynamic framework is restricted to the interiors. Hence, the
zero-gradient condition is utilized to hold back sharp nonhydrostatic
pressure gradients at the open boundaries. With the above boundary
conditions, this linear system (Eq. 28) can be solved via the Krylov subspace
method with PETSc's assistance on parallel computers under the standard
MPI-based framework (Balay et al., 2020). Also, a highly efficient method
needs to be devised to precondition the huge and sparse
matrix

**A**. Here, the multigrid preconditioner (Smith et al., 1996) and flexible generalized minimal residual algorithm (Saad, 1993) are employed in numerical validation experiments in this paper to minimize computational costs.

In this section, we present a series of ideal numerical validation experiments to explore the correctness and compatibility of nonhydrostatic algorithms together with ORCTM. In allusion to the internal solitary wave dynamics, these test cases range from laboratory-scale cases in an enclosed tank to field-scale ones like the northern South China Sea with open boundaries. The first case is the lock-exchange problem as the preliminary validation. The second to fourth cases are designed to explore the nonlinear evolution of internal solitary waves induced by their interactions with the changing terrain. The last one is the generated nonlinear internal wave case in a double-ridge environment analogous to the Luzon Strait, which aims at the generation and disintegration of nonlinear internal waves to examine the effectivity of the open tidal forcing condition module under the nonhydrostatic algorithms. Analyses of all test experiments above indicate that ORCTM can reproduce nonlinear and nonhydrostatic internal solitary waves in different oceanic environments, which shows the robustness and reliability of this nonhydrostatic ocean model.

## 3.1 The lock-exchange problem

When shear currents flow between two different density fluids, the
Kelvin–Helmholtz instability (hereafter K-H instability) will appear to
cause turbulent diapycnal mixing (Lawrence et al., 1991; Cushman-Roisin,
2005). The perturbation on the interface gradually develops and stimulates
numerous small eddies due to energy dissipation. The magnitude order of
vertical flow is comparable to the horizontal one, so the nonhydrostatic
effect matters throughout the whole process. We set a rectangular enclosed
tank separated by a vertical board in the middle at the *x*-axis origin. Both
sides of the tank are separately filled with two different density fluids in
Fig. 1a. The gravitational adjustment will proceed when the central board is
disengaged just like a lock gate. Here, we refer to the previous
configurations (Härtel et al., 2000; Fringer et al., 2006; Lai et al.,
2010) as a 2-D problem. The horizontal length *L* is set to 50 cm, and the
static water height is 10 cm without topographic change in the tank. The
grid resolution is 0.001 m in the horizontal and vertical directions. Several
sensitivity experiments were explored to reduce the dissipations out of
solid boundary friction, so the bottom friction coefficients are finally set
to zero; both *A*_{V0} and *D*_{V0} in Eqs. (11) and (12) are $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ m^{2} s^{−1}. In addition, water density
averages are calculated based on the prescribed salinity difference on the
left and right sides of the tank: *ρ*_{l}=1023.05 kg m^{−3} and *ρ*_{r}=1026.95 kg m^{−3}.

The K-H instability process grows rapidly with good eddy reconstruction
and outstanding wave breaking. In contrast, in that model configuration, we
also run the same configuration experiment above but under the hydrostatic
balance scheme. Figure 1b and c show the results of density *σ*
(define $\mathit{\sigma}=\mathit{\rho}-\mathrm{1000}$ kg m^{−3}) at the same time under the
hydrostatic and nonhydrostatic balance framework. The comparison proves
that the K-H instability cannot proceed resulting from the inapplicability
of the hydrostatic balance. The perturbation on the density interface is so
tiny that the density fronts cannot evolve in the upper and lower layer, so
the mixing caused by the overturning and shear is too weak to be seen. On
the contrary, via the nonhydrostatic scheme, the eddies can proliferate with
energy dissipation due to the associated shear on the perturbation,
vigorously mixing the high- and low-density water on the interface. More
specifically, the energy is transmitted to the small-scale eddies across the
density fronts due to dispersion and nonlinearity.

The evolution process of K-H instability is shown in Fig. 2. It is out of gravitational adjustment that the density front movement is accompanied by heavy water in the bottom and light water in the upper level moving to the left and right, respectively, causing a velocity shear field and clockwise rotating interface in Fig. 2a. The shear strength gradually increases until breaking the critical point of restoring force that depends on the density gradient, and later a series of eddies grows from the middle to both sides of the tank with the turbulent rolling and overturning. These eddies mix the water body with high density at the bottom and the upper one with low density, forming multiple considerable density mixing areas in Fig. 2b and c. When the bottom density flow is reflected on the left wall, a similar adjustment process begins to develop in reverse of Fig. 2d and e, but the strength of subsequent eddies is significantly weakened due to energy dissipation.

These density distributions display the generation of density fronts and numerous eddies throughout the gravitational adjustment process. Based on the point of energy dynamics, the gravitational potential energy (PE) is converted to kinematic energy (KE) for the water parcel, while the total energy dissipates continuously in the tank. Here, KE and PE of the entire 2-D tank are calculated from the following formulas.

The three curves show the fluctuation of PE, KE, and total energy during the K-H instability simulation in Fig. 3. The PE and KE correspond to the maximum and zero due to the initial density distribution and static field in Fig. 3a. Afterward, the PE declines sharply with an opposite change in KE. Both rates of change are almost the same based on the curve slopes, which demonstrate that PE is converted to KE, reaching mutual peaks of about 9.5 s at the end of the first gravitational adjustment. From then on, both of them still maintain the opposite trends with an oscillation of roughly 25 s. It is worth noting that all kinds of energy exhibit a downward trend, with their oscillation period increasing steadily due to energy dissipation so that KE will drop to zero, and PE and total energy (PE + KE) will reach constant values in the end. The results above are equivalent to previous works (Härtel et al., 2000; Fringer et al., 2006; Lai et al., 2010), implying the correctness of the nonhydrostatic dynamic module.

## 3.2 Internal solitary wave in a tank

Internal solitary wave activities are ubiquitous in the ocean, with strong
nonlinearity and nonhydrostatic effects. Laboratory experiments are usually
carried out to study the ISWs to make up for the deficiencies of field
observations. Numerical ISW experiments in a laboratory-scale background
need to be combined simultaneously (Grue et al., 2000). We follow the
previous experimental configuration (Ma et al., 2020). A schematic diagram
of the ISW experiment is given in Fig. 4. The tank length is 2.0 m with the
*x*-axis origin located on the left; the static height is 10 cm without
topographic change; the horizontal and vertical resolutions are $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ m; both bottom friction and drag
coefficients are set to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$ with the effect of a fairly
robust friction to the ISW; *A*_{V0} and *D*_{V0} are the same as in the experiment configuration in Sect. 3.1. Here, a gravity
collapse method is used to generate the depression ISW. Specifically, the
low- and high-density fluids initially fill the upper and lower layers of
the tank with the collapse area on the left side. The collapse height and
width are 5.0 and 4.0 cm. Water density averages are calculated in the
upper and lower layer with *ρ*_{1}=1003.62 kg m^{−3} and *ρ*_{2}=1026.95 kg m^{−3}. Additionally, a diagnostic module is employed to
characterize the high-frequency variation. The high-frequency outputs are
positioned at points *x*=0.4, 0.8, 1.2, and 1.6 m with a time interval of
0.05 s.

Figure 5 distinctly illustrates the evolution of the ISW packet in the tank
based on the pycnocline fluctuation. The isopycnic of 1026 kg m^{−3} can
characterize the maximum strength of the depression ISW in Fig. 5a. The eastward
starting wave packet originating from the west gravity collapse area
comprises the depression heading wave and several tail waves whose
amplitudes decrease successively. The heading wave with the maximum
amplitude propagates much faster than the tails behind so that the distance
expands promptly between them. As shown in Table 1 about the heading
wave characteristics at the four locations, we find the wave amplitude with
almost little change and then a slight fluctuation but no more than 0.1 cm
after *x*=0.8 m. The quantitative evaluation of the wave speed based on the
slope of the blue dashed area in Fig. 5b reveals that the wave speed
increases slowly after *x*=0.2 m but with its increment less than 0.01 m s^{−1}. This indicates that the starting ISW packets are still at the
stage of gravity adjustment before arriving at *x*=0.2 m and then
propagating to the east steadily in our simulation. Also, the
characteristic westward-reflected waves (the blue line in Fig. 5a) with the
larger amplitude prove that wave–wave interactions happen between the
reflected and starting tail waves.

We select a snapshot result for characteristic verification shown in Fig. 6
when the heading wave arrives at around *x*=0.8 m. The strongest horizontal
velocity of the depression wave is 0.023 m s^{−1}, and the vertical flow
can reach up to 0.0065 m s^{−1}. The characteristic velocity fields are in
line with the clockwise structure of a theoretical depression internal
solitary wave. Furthermore, the nonlinearity $\mathit{\epsilon}=a/h$ and dispersion $\mathit{\mu}={\left(h/\mathit{\lambda}\right)}^{\mathrm{2}}$ are calculated at the different
locations in Table 1, where *a*, *h*, and *λ* are the amplitude, water
height, and characteristic wavelength. The KdV (Korteweg–de Vries) model (Benjamin, 1966)
described in Appendix B is utilized to predict theoretical waveforms at the
four locations. The comparison depicted in Fig. 7 demonstrates that the
results are more consistent with the KdV model than the m-KdV (modified KdV) model.
According to the nonlinearity *ε* from Michallet and
Barthélemy (1998), small- and large-amplitude ISWs can be classified
when *ε*<0.05 and *ε*>0.05,
respectively, whereas the application of the KdV model requires a balance
between weak nonlinearities and dispersion (Ono, 1975), which namely
needs to satisfy this condition: $\mathit{\mu}=O\left(\mathit{\epsilon}\right)\ll \mathrm{1}$. Despite the large-amplitude waves simulated from our
model with *ε*>0.05, the nonlinearity and dispersion
are of the same order and small enough that the heading wave can be deemed
under weak nonlinearity. This can explain why the waveforms are
better described by the KdV model. Therefore, analyses of the theoretical
model indicate that the simulation of internal solitary waves can be
fulfilled authentically using our nonhydrostatic model.

## 3.3 Internal solitary wave shoaling on a Gaussian terrain

Based on the experiment configuration in Sect. 3.2 (also called Exp. 3.2), a slowly varying terrain is implemented to explore the nonlinear
evolution of internal solitary waves in Sect. 3.3 (also called Exp. 3.3), especially the wave shoaling. As shown in Fig. 8, the left half of the
Gaussian curve is reserved as the slope-shelf terrain starting between *x*=1.0 and 1.3 m with a height of 5.0 cm, and then the water depth remains
unchanged from *x*=1.3 to 2.0 m corresponding to the shallow-water zone.
High-frequency outputs are acquired during the climbing process of ISWs
at points *x*=0.4, 0.8, 1.0, 1.2, 1.3, 1.4, 1.6, and 1.8 m with the same
output interval as Exp. 3.2.

The evolution of internal solitary waves with varying topography is
displayed in Fig. 9. The heading ISW holds a stable packet at *x*=0.4 m and
initiates shoaling after reaching *x*=1.0 m. Afterward, the heading ISW
undergoes topographic change so that the speed of the wave trough is
less than the wave rear. Consequently, the contrasting effects on the wave
front and wave rear contribute to the former gentle sloping but the latter
gradual steepening, which shows a similarity with Vlasenko et al. (2002).
Then the closed isopycnic contour mirrors the backward overturning and
rolling due to the wave breaking at *x*=1.2 m in Fig. 9a. Apart from the
wave breaking process above, it is also found in Fig. 9b that the reflected
waves propagate back to the deep-water zone from *x*=1.2 m. In other words,
both wave breaking and refection substantially attenuate the original
depression ISW energy. When arriving east of *x*=1.2 m, the original
depression wave is past the critical point where the upper layer is thicker
than the lower one in Fig. 10a, so an elevation wave springs up in the wave
rear. The elevation wave then continues to propagate eastward, which leads
to increasing accumulation of high-density water in the upper water in
the right region close to the wall of the tank. Hence, a new collapse area
between *x*=1.8 m and the east wall comes into being where the thickness of
the upper layer is larger than the lower layer. Ultimately, the westward-reflected waves, including a series of elevation tail waves, are released at
*x*=1.6 m. In detail, the first elevation is the leading one with a
rank-ordered structure in the rear. After reaching the deep-water zone left
to 1.3 m, the wave rear begins to steepen and sink, and a depression wave
forges behind the wave rear. Namely, the soliton wave passes the critical
point inversely due to the wave deepening.

For further exploration of the evolution of the depression wave, the
distributions of the vorticity $(\mathit{\zeta}=\partial w/\partial x-\partial u/\partial z)$ with velocity vector
are depicted in Fig. 10. The depression wave core features negative vorticity
with an anticyclonic velocity structure before reaching the shelf
topography. When the ISW approaches the top of the slope in Fig. 10a and b, the vertical shear increases promptly and strengthens the positive
vorticity at the bottom. Then the backward overturning springs up between
*x*=1.2 and 1.3 m, marking the ISW entering the breaking instability stage
(Helfrich and Melville, 1986) due to the shoaling. At this time, even though
wave breaking and reflection render partial wave energy dissipation,
a fraction of the depression wave can reach the shallow-water zone,
leaving a cyclonic vortex behind above the slope shelf in Fig. 10c. This
partial soliton wave is adjusted instantaneously when the upper layer
thickness is more significant than the lower in light of the boundary of
the negative vorticity area in Fig. 10d. As a result, after the reverse
situation occurs, the elevation wave begins to emerge at the back of the
original wave. Its core corresponds to the positive vorticity with a
cyclonic velocity structure. In addition, the vortex from the wave breaking
weakens slowly and motivates numerous small-scale waves with high wavenumber
propagating to both sides in Fig. 10e and f, which is consistent with the
propagation characteristics of the reflected waves near *x*=1.2 m in Fig. 9b.

It is also worth highlighting the evolution of the reflected westward
waves. We also visualize the process of the second reverse situation due to
the wave deepening in Fig. 11. It can be noticed that there is a leading
elevation wave at *x*=1.4 m followed by a series of rank-order waves
exhibiting a likewise sinusoidal variation. They propagate together to the
deep-water zone with the wave crest corresponding to positive vorticity.
Particularly, the wave train is considered approximately linear based on
the alternating positive and negative vorticity regions, since the cores of
these waves are located almost in the middle layer where the nonlinear
parameter *α* is close to zero in terms of the KdV model. As the water
depth becomes deeper, the crest of the elevation wave gradually grows down
and flattens with the wave rear sinking. The original elevation cannot be
maintained in the deep water, transforming into a depression wave with the
velocity fields adjusted accordingly.

The ISW is not stable enough to coincide with the KdV model after passing
the critical point into the adjustment stage. Hence, we select the two types
of soliton results for verification before the reverse situation occurs. The
comparison results between theoretical and numerical model are illustrated
in Fig. 12 at *x*=0.8 and 1.4 m before the wave shoaling and deepening,
respectively. We can find that the depression waveform conforms to the KdV
model results before climbing the slope, whereas the elevation is closer to
the m-KdV model. Compared with *ε*=0.233 at *x*=0.8 m,
the interaction between the ISW and the shoaling topography renders a
stronger nonlinearity (*ε*=0.331) of the elevation
heading wave in the shallow water. Namely, the larger wave amplitude ratio
in the shallow-water results can be characterized with m-KdV theory, which
compares well with the conclusions of Michallet and Barthélemy (1998)
in a satisfactory way.

## 3.4 Internal solitary wave breaking on a slope.

To further characterize a complete breaking and dissipative process of ISWs,
we set a linear slope identical to Michallet and Ivey (1999). As is shown in
Fig. 13, the tank length is 2.0 m; The height is 15 cm with the linear
terrain placed on the east side. The model configuration (i.e., spatial
resolution and viscous coefficients) is identical to Exp. 3.2, which can
ensure the same time step according to the Courant–Friedrichs–Lewy (CFL)
condition, and the depression ISW is about to be dissipated due to
increasing bottom friction at the shelf break. In contrast with Bourgault
and Kelley (2004), water density averages are calculated to be *ρ*_{1}=1000.01 kg m^{−3} and *ρ*_{2}=1047.00 kg m^{−3} in the upper and lower layers. Via several sensitivity experiments
about collapse area, the amplitude of the depression wave can reach
approximately 2.8 cm when the collapse height is 9.0 cm with its width
identical to Exp. 3.2. Although the stimulated wave strength is slightly
greater than the results from Bourgault and Kelley (2004) due to the
different wave generation methods, it is predictable that the breaking of
the larger-amplitude ISW will be more dramatic with a prominent performance
for model verification.

The associated density and velocity fields produced by the depression ISW
are presented in Fig. 14 at *t*=15 s before wave shoaling. The horizontal
velocity is about 3.0 cm s^{−1} at the surface and varies up to 3.5 cm s^{−1} at the wave core. Meanwhile, the vertical velocity distribution
presents a double-core structure reaching ±0.8 cm s^{−1}. The
unique anticyclonic velocity characteristic just like an eastward rolling
wheel is consistent with the model results of Bourgault and Kelley (2004).
We select the four moments of the evolution of wave shoaling illustrated in
Fig. 15. In addition to wave breaking accompanied by the waveform steepening
in the rear, a significant density front rolling in the wave front evolves
along the linear slope during the overall shoaling process in Fig. 15a and b. Specifically, while the depression wave continues to get closer to the
shallow zone, the effect of bottom friction can maintain the vertical shear
and increase the potential energy, which intensifies the diapycnal mixing
and dissipation on the density interface. Then the wave-induced diapycnal
flow contributes to high-density water under the interface transported
continuously to the shallow zone in Fig. 15c. On the other hand, there is
another pronounced peculiarity in Fig. 15d compared to Exp. 3.2. A few
small-scale eddies emerge along with the sheared interface due to the shear
instability.

To further evaluate and validate the wave breaking process, we compare the velocity field distributions with the observation results via PIV (particle image velocimetry) technology from Michallet and Ivey (1999) and nonhydrostatic numerical experiments from Bourgault and Kelley (2004) in Fig. 16. Accordingly, when the depression wave arrives over the slope, its depression waveform and anticyclonic flow field are modulated by the topographic shoaling to flatten the wave front and enhance the downward current along the slope due to the bottom friction. Meanwhile, a smaller cyclonic eddy appears and clings to the slope under the steepened wave rear in Fig. 16a. As the deformed depression wave persists in shoaling, the cyclonic eddy is reinforced and extends its scope of influence, resulting in a strong overturning from near the bottom layer to promote the wave steepening in Fig. 16b, which presents good agreement with the results from Exp. 3.3. Afterward, the anticyclonic flow structure has been ruined since bottom friction commences, hindering the current down the slope. In contrast, the coverage of the cyclonic eddy continues to expand and moves the shallow zone with the waveform distorted further. All the above nonlinear processes are similar to previous laboratory and model results. Our nonhydrostatic model can also resolve the nonlinear evolution of the internal solitary waves at shelf break with high enough accuracy.

## 3.5 Nonlinear internal waves in a double-ridge system

The last validation experiment is to examine the generated nonlinear internal waves via tidal flow over the varying topography. We set up an underwater double-ridge system comparable to the Luzon Strait in the northern South China Sea (SCS), where the largest internal solitary waves in the world can exist (Huang et al., 2016). This validation case is a 2-D problem for the reduction of computational resources as well. The topography in this double-ridge system is fitted approximately with the Gaussian function given as

In Eq. (32), *H*(x) is the water depth; the height of the east ridge
and west ridge (*h*_{e} and *h*_{w}) is 2500 and 1300 m in sequence with an
interval and widths of 100 km, which is similar to the fundamental
topographic characteristics in the Luzon Strait. As shown in Fig. 17a, the
static water height is 3000 m; the east ridge and west ridge (hereafter ER and WR)
are located at the coordinate origin and $x=-\mathrm{100}$ km; the
horizontal and vertical grid resolutions are uniformly 200 and 10 m
throughout; *A*_{V0} and *D*_{V0} in Eqs. (11) and (12) are set to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ m^{2} s^{−1}; the bottom friction coefficients are both
$\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$. As for the tidal categories, the generation
of semidiurnal ITs and the modulation effect of diurnal ITs in the Luzon
Strait determine the evolution of the larger-amplitude ISW packet in the
northern South China Sea (Buijsman et al., 2010a; Zeng et al., 2019), so we
define the *M*_{2} and *K*_{1} tidal current amplitudes as 5.0 and 4.0 cm s^{−1} corresponding to the semidiurnal and diurnal components at the open
boundaries; the sponge thicknesses *L*_{sp} of the west and east
boundaries are both approximately 40 km, and the timescale coefficient *τ* is set to 500 s. These model configurations in the validation
experiment are analogous to the control test from Li (2014) and Zhang et al. (2011) to reproduce the major structures of NIWs in the South China Sea.
In addition, to simplify the background environment, we also use horizontally
uniform stratification as the initial field for our model. Here, the
representative stratification in Fig. 17b to d stems from the GLORYS12V1
reanalysis product in CMEMS (Copernicus Marine Environment Monitoring
Service). The initial field is based on the spatial mean around the source
of generated ISWs in the Luzon Strait during the summer of 2011, since
large-amplitude ISWs are observed during this period on the SCS continental
shelf (Ramp et al., 2019) and the strong thermocline structure in summer is
conducive to the formation of baroclinic tides in the Luzon Strait (Zheng et
al., 2007; Buijsman et al., 2010b;). Additionally, the slope criticality *γ* (Gilbert and Garrett, 1989; Shaw et al., 2009) no less than 1 is
usually essential with the formation of linear internal waves:

in which *ω* is the tidal angular frequency, *N*^{2} is
the buoyancy frequency squared, and the Coriolis parameter *f* is set to zero
for the Earth's rotation, which is neglected due to the 2-D environment. Around the east
ridge *γ* is always larger than unity regardless of the *M*_{2} and
*K*_{1} tide, which means the east ridge belongs to the supercritical
topography. Therefore, it is predictable to generate internal waves due
to the interactions with barotropic flow over the east ridge. We run the
model for 10 d from an initial static field. The diagnostic module is
also used to characterize the high-frequency variation with the output
interval of 1 min at $x=-\mathrm{250}$, −350 km.

Figure 18 shows the maps of horizontal baroclinic velocity
$u{}^{\prime}=u-U$, where *u* is the total velocity and
*U* is the barotropic flow velocity. From the characteristics of the source
field, the generation of an internal tide beam propagating
eastward and westward centered from the eastern side of the east ridge is found. The
eastward barotropic flooding current flows continuously over the east ridge
with a maximum barotropic current up to 0.0531 m s^{−1}. A significant
hydraulic jump can appear with the isotherm fluctuation up to roughly 200 m
on the eastern side, which indicates the formation of lee waves to a certain
extent. The above internal wave generation due to tide–topography interactions
can be described with the below nondimensional parameters at the source: (1) the tidal excursion parameter $\mathit{\epsilon}={U}_{\mathrm{0}}/L\mathit{\omega}$, which can be
associated with the generation of an internal tide beam under critical or
supercritical topography, where *U*_{0} is barotropic current
amplitude from the far field and *L* is the characteristic length for
topography (Garret and Kunze, 2007; Chen et al., 2017); (2) the Froude
number $Fr=U/c$ and its topographic form
$F{r}_{z}=\mathit{\omega}/N(\mathrm{d}H/\mathrm{d}x)$, in which *c* is the mode-1 linear speed for the eigenvalue
problem (Legg and Adcroft, 2003; see Appendix B). Specifically, Legg and
Klymak (2008) found that the nonlinear hydraulic jump will develop with lee
wave generation when $F{r}_{z}<\mathrm{1}/\mathrm{3}$. It is worth noticing that the
tidal excursion far less than unity agrees with the formation of the linear
internal tide beam on the critical or supercritical topography but cannot
ensure the formation of lee waves altogether. For instance, the lee
waves remain strong in the Luzon Strait despite the tidal excursion under
unity (*ε*≈0.4) in previous model results (Buijsman et
al., 2010b). The tidal excursion parameter *ε* and the Froude
number *F**r* are estimated to be 0.025 and 0.018 in Fig. 18. That demonstrates
that the multi-modal baroclinic tides and upstream propagation of internal
waves will generate around the source field when the subcritical barotropic
current flows over the east ridge. Furthermore, the maximum topographic
Froude number is just 0.3362 around the east ridge with the approach to the
regime transition value $\mathrm{1}/\mathrm{3}$, which ensures that the nonlinear hydraulic jump
can grow with lee waves on the east of the east ridge. All of the above can
explain the generation of the internal tide beam and hydraulic jump well in
our simulation and confirm the mixed tidal lee wave regime in the Luzon
Strait (Chen et al., 2017).

After the westward internal tide beam emitting from the east ridge reaches the sea surface and reflects into the deep sea, the partial downward internal tide beam can propagate to the top of west ridge below 1500 m depth and reflect into the upper layer again. Between the double ridges, such a significant portion of beam energy captured by the pycnocline waveguide together with the upstream influence can strengthen the westward-propagating internal wave energy in Fig. 18b, which can trace back to the source of the internal solitary wave packets in the far field. However, the strong dissipation for the high modal internal waves contributes to the vanishing of the internal tide beam structure and allows the nonlinear evolution of low-mode baroclinic tides. The low modal internal solitary wave packets can grow and propagate westward from $x=-\mathrm{150}$ km, marking the disintegration of the multi-modal nonlinear internal wave energy. Specifically, the first-mode ISW packet emerges from $x=-\mathrm{250}$ to −200 km. Meanwhile, the second-mode ISW between $x=-\mathrm{350}$ and −300 km performs the convex wave packet.

We can acquire the propagation characteristics of these ISWs via analyzing the global temperature time series at 400 m depth. As is illustrated in Fig. 19, the second-mode ISW packet propagates slower, and its strength is much weaker than the first-mode wave one. Also, it can be determined that the two first-mode wave packets can propagate westward in 1 d, one of which is stronger with the structure of several tail waves, and the other is almost solitary and weak. These two types of first-mode wave packets refer to type-a and type-b waves (hereafter a-wave and b-wave), respectively, in the northern South China Sea (Ramp et al., 2004). In addition, their occurrence time can be connected to the ebb of the eastward flood current around the east ridge. These simulated results in the strength and timing prove that a- and b-waves originate from the double ridge in Luzon Strait (Ramp et al., 2004, 2019; Zhao and Alford, 2006). Additionally, the relatively weak second-mode concave wave can be found distinctly following the a-wave from the west of −300 km. To sum up, the multi-modal baroclinic tide structures from the double-ridge system can propagate to the far fields. The low-mode internal waves gradually perform the corresponding ISWs due to nonlinear enhancement, which displays good agreement with the other two-dimensional experimental results (Buijsman et al., 2010a, b; Vlasenko et al., 2010).

To evaluate the comparison between the numerical ISWs with internal wave
theory, we select the results of the first-mode ISW at $x=-\mathrm{250}$ km. In Fig. 20, it is found that a first-mode ISW packet including three
tail waves arrives at the position after 10:00 on the sixth day. The maximum
fluctuation of the first-mode ISW packet can reach 206 m located between 650
and 900 m water depths. The westward horizontal baroclinic velocity
associated with the wave packet prevails above 200 m with a maximum strength
of roughly 1.41 m s^{−1}, and the corresponding downwelling region is
located between 200 and 1500 m depths with the strongest downward velocity
up to 0.22 m s^{−1}. According to the sea surface height gradient (SSHG,
SSHG is defined $\sqrt{(\mathrm{\nabla}\mathit{\varsigma}{)}^{\mathrm{2}}}$) in
Fig. 20d, the average propagation speed of this wave packet is approximately
3.17 m s^{−1} based on the slope of the SSHG contour. Moreover, we solve the
Taylor–Goldstein equation (Miles, 1961; Liu, 2010; see Appendix B) 10 min before this wave packet reaches $x=-\mathrm{250}$ km, and the normal mode
of vertical velocity is subject to the rigid-lib boundary condition. We
found that the location of the maximum modal function is 710 m, in agreement
with the model results in Fig. 20e. However, the propagation speed is
greater than the first-mode linear result of 2.69 m s^{−1}, which is
probably attributed to the underestimated effect in linear theory.
Therefore, the KdV model is also utilized to analyze the depression wave
packet. The nonlinear and dispersion parameters are $-\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{3}}$ s^{−1} and 2.4×10^{5} m^{3} s^{−1}, which indicates that
the theoretical depression wave is consistent with the simulated results
(Helfrich and Melville, 1986). Nevertheless, the theoretical nonlinear
velocity of about 2.88 m s^{−1} is slightly lower than the simulated
results. It is probable that the increasing nonlinearity with the steepening
of internal tides ultimately leads to the larger propagation speed of this
first-mode ISW packet.

It is also noticeable that the multi-modal internal solitary wave field is
generated and strengthened gradually due to nonlinear enhancement. In
Fig. 21a, we can recognize the distinct ISW packets from the isotherm
displacement that refers to type-a waves, second-mode waves, and type-b waves from
the source to the far field. The a-wave packet features the most substantial
strength with tail waves when its vertical excursion induced by the heading
wave can reach up to 120 m. In contrast, the weaker b-wave contains one
depression soliton in the west to $x=-\mathrm{400}$ km. They both originate from a
multi-modal internal tide caused by the tide–topography interactions in the
double-ridge system, but the b-wave is more associated with the west ridge
(Buijsman et al., 2010a; Zeng et al., 2019). Between a- and b-waves, there is
a second-mode ISW packet evidently classified as having the structure of a concave wave
whose upper and lower isotherms fluctuate downward and upward. The maximum
isotherm fluctuations are located in roughly 180 and 1000 m depths and
can reach up to −57.2 and 140.6 m. The propagation speed of this second-mode
ISW is about 1.36 m s^{−1} from the SSHG slope in Fig. 20d. It is
predictable that the a-wave packet will follow the second-mode signal due to
the more considerable speed. Figure 21b and c show the second-mode ISW
packet and related velocity field time series at $x=-\mathrm{350}$ km. The
horizontal baroclinic velocity field has a sandwich-shaped vertical
structure, and the maximum of 0.42 m s^{−1} is located in the middle layer
between 200 and 600 m. The strength of baroclinic velocity with a small
average of 0.2 m s^{−1} is distinct from the stronger first-mode ISW
packet above 200 m. Additionally, a double-peak structure performs in the
vertical velocity field, and it is distributed at the depths of 150 and 1000 m where the strength in the deep layer is stronger than the upper, resulting
in a minor isotherm fluctuation above 200 m. Here, the Taylor–Goldstein
equation is also solved to acquire the eigenfunction of the vertical
velocity. In Fig. 21d, the second-mode eigenvalues have two vertical peaks
whose depths correspond to 150 and 1070 m, with the latter strength stronger
than the former, and the corresponding phase speed is about 1.34 m s^{−1}.
In summary, first- and second-mode internal solitary waves as the leading
carriers can transfer baroclinic tidal energy from the source to far
fields until dissipating thoroughly. The multi-modal solitary wave field
conforms with the previous two-ridge experimental result using MITgcm
(Vlasenko et al., 2010). Internal wave theoretical models can compare
well with the distribution of stimulated results in our nonhydrostatic ocean
model, demonstrating an overall good performance of characterizing the
nonlinear evolution of multi-modal baroclinic tides.

The main focus of this paper is to introduce a new oceanic regional nonhydrostatic circulation and tide model (ORCTM) which is rooted from the MPI-OM and aims to characterize the internal solitary wave processes of real oceans, such as in the northern South China Sea. We developed and implemented the nonhydrostatic dynamics and open boundary module under the original global hydrostatic framework of MPI-OM. Based on the fractional step and finite-difference methods, ORCTM involves the three-dimensional fully nonlinear momentum equations under the Boussinesq fluid. The three-dimensional Poisson equation subject to different boundary conditions needs to be solved before the pressure correction method is employed to acquire the velocity field corrected via nonhydrostatic pressure gradient force. In order to match the nonhydrostatic algorithm and realize larger-amplitude ISW simulations in an ocean-scale case, an exponential relaxation term is implemented to the control equations through the sponge layers as the open boundary condition.

A series of two-dimensional ideal numerical experiments associated with the nonlinear evolution of the internal solitary waves and baroclinic tides is devised to verify this nonhydrostatic ocean model. Here, the results of the validation experiments are in accord with the theoretical framework of the nonhydrostatic dynamics and demonstrate that ORCTM can successfully characterize the generation, propagation, and dissipation of internal solitary waves in laboratory-scale cases. Specifically, the reverse situation due to wave shoaling and deepening can be depicted completely when considering the topographic change. Meanwhile, the stimulated internal solitary wave conforms with the previous numerical experiment and the direct observations in the laboratory test. Also, ORCTM can capture the density fronts with the cyclonic eddy induced by the wave breaking, which shows good stability and high enough accuracy. Furthermore, based on the real topographic features in the Luzon Strait of the northern South China Sea, analyses of the validation experiment indicate the multi-modal structure of baroclinic tides in the double-ridge system. The nonhydrostatic ocean model ORCTM is proven to be able to reproduce the life cycle of multi-modal ISWs induced by tide–topography interactions in the Luzon Strait and precisely capture the alternation process of type-a and type-b internal solitary wave packets. The first two-mode ISW structure compares well with the internal wave theoretical model.

Even though these validation experiments have a strong resemblance to other
nonhydrostatic model results (Bourgault and Kelley, 2004; Berntsen et al.,
2006; Lai et al., 2010), some distinctions in grid structure or numerical
methods may have an opposite impact, especially when predicting a particular
nonhydrostatic dynamics process. Berntsen et al. (2006) indicated some noisy
structures near the bottom layer due to numerical errors of finite-volume
treatment when predicting internal solitary wave breaking via MITgcm
(Marshall et al., 1997a, b). They found that the Bergen Ocean Model (BOM) can avoid this
problem with a sigma coordinate, whereas MITgcm needs a high-order filter to
suppress the noise. However, artificial flow usually emerges and has a
negative influence on ISW breaking simulations due to internal
baroclinic pressure errors in the sigma coordinate. These require model
users to have a refined grid when encountering an area of changing topography.
Compared to the nonhydrostatic unstructured grid finite-volume coastal ocean model (FVCOM) (Lai et al., 2010) and BOM (Berntsen et
al., 2006), ORCTM is based on the finite-difference method and owns a
*Z* coordinate, which has the capability to avoid the above errors. These
numerical methods and validation experiments demonstrate that ORCTM is able
to approach or reach an acceptable or better level for a nonhydrostatic ocean
model for ISW simulation.

The simulation of internal solitary waves can mirror the macroscopic structure and assist with the implementation of in situ observations. It is noticed that the predictability of nonlinear internal wave characteristics relies on the model performance and external conditions such as realistic stratification, bathymetry, and background circulation. Another advantage of ORCTM is the usage of an orthogonal curvilinear mesh grid in the horizontal direction. It is competent enough to maintain the small-scale nonhydrostatic dynamics well-resolved in the concerned region via mesh refinement. Particularly, constructing practical and reliable background fields via a nested technique remains the way forward for ISW simulations in the oceanic environment. Enhancing the fidelity of ISW simulation remains challenging. Nevertheless, it can be concluded that our regional nonhydrostatic ocean model is a good choice for oceanography scientists interested in internal wave research and numerical prediction.

According to the idea of fractional steps (Chorin, 1968; Press et al., 1988), a pressure correction method on the nonhydrostatic dynamic component is employed to calculate the intermediate velocity over the original hydrostatic balance scheme (Fringer et al., 2006; Lai et al., 2010). If the flow is close to the hydrostatic balance, the pressure of the nonhydrostatic part will be so slight that the correction plays a minor role. The key to the nonhydrostatic dynamics module is to solve the Poisson equation below.

The right-hand side (RHS) of Eq. (A.1) is the divergence about the intermediate velocity as a source or sink term. Here, based on the definition of divergence, the three components calculated directly at each cell are specified in the three orthogonal coordinates as follows.

Here, $i,j,$ and *k* are the indices of increasing eastward, northward, and
downward along the *x*, *y*, and *z* axis, respectively; *z*=0 is defined on the
undisturbed sea surface by means of local Cartesian coordinates;
$\stackrel{\mathrm{\u0303}}{u}$, $\stackrel{\mathrm{\u0303}}{v}$, and $\stackrel{\mathrm{\u0303}}{w}$ are the intermediate velocity;
Au, Av, and Aw represent the six faces area of a cell in the $i,j,$ and *k*
directions; and Ω is the volume of a cell. These grid descriptors are
defined as follows.

DX, DY, and DZ represent the spacing difference between the
adjacent grid cells in the *x*, *y*, and *z* axis. The suffixes associate *u*, *v*, and
*w* at the cell face center and *p*^{′} at the body center. Compared to the finite-difference method, the definition of the divergence of a cell is more
accurate and reliable, especially when adjacent to the solid boundaries for
the RHS calculation. The left-hand side (LHS) of this equation is
discretized horizontally on the Arakawa-C grid (Arakawa and Lamb, 1977)
using the central difference method with second-order accuracy. The
vertical discretization is the same as the Max Planck Institute Ocean Model
(Marsland et al., 2003), wherein the bottom grid has the capacity of the
partial filled cell to adjust the vertical distance for fitting into the
realistic terrain (Marshall et al., 1997b). We can acquire the following
finite-discrete equation about seven cells for nonhydrostatic pressure
perturbation as

where the coefficients of the discretized LHS are given as follows.

Invoking the boundary conditions (Eq. 29) and Eqs. (A.6) to (A.7), the discretized Poisson equation with seven cells can be derived with the matrix form below:

where **A** is a sparse and definite-positive matrix with seven
diagonals; ${{p}^{\prime}}_{\mathrm{nh}}$ and ** B** are the column vectors
with a size of all cell numbers $Nxyz=Nx\times Ny\times Nz$ in the model domain, where

*N*

*x*,

*N*

*y*, and

*N*

*z*are the cell numbers in the $i,j,$ and

*k*directions. Actually, the sparse matrix

**A**cannot easily be handled directly with a size of

*N*

*x*

*y*

*z*×

*N*

*x*

*y*

*z*, which hence needs to be designed with greater efficiency as a precondition. To apply the nonhydrostatic model to the real oceanic environment on the original model base, the Portable, Extensible Toolkit for Scientific Computation (PETSc) library is implemented into the nonhydrostatic dynamic module. We apply the numerical Krylov subspace methods for the matrix solvers under an MPI-based framework (Balay et al., 2020). Here, the flexible generalized minimal residual (FGMRES) method (Saad, 1993) is applied to solve this problem in conjunction with a multigrid preconditioner (Smith et al., 1996) for the sparse matrix before iteration. Thus, the nonhydrostatic pressure can be computed with these methods. It should be emphasized that the nonhydrostatic and hydrostatic dynamics modules remain independent of each other and not contradictory. The nonhydrostatic dynamics module will make up for the deficiency of the hydrostatic module only considered in this model, which means the nonhydrostatic and hydrostatic simulations can be simultaneous in this model. In other words, the nonhydrostatic dynamics can be fulfilled economically in harmony with the original numerical framework.

Based on the shallow-water approximation, a small-amplitude internal solitary wave whose amplitude compared with the total depth is small enough can be described by the classical two-dimensional Korteweg–de Vries (KdV) equation given as follows (Apel et al., 2007).

Considering a two-fluid stratification system is more appropriate for the
experiments in Sect. 3.1–3.3. *ρ*_{1} and *ρ*_{2} are the upper
and lower densities corresponding to the thickness *h*_{1} and
*h*_{2}; *x* is the horizontal coordinate. Several parameters can be
written here as (Benjamin, 1966; Wessels and Hutter, 1996)

where nonlinear and dispersion parameters (*α* and *β*,
respectively) can represent the soliton polarity; *c* is the linear velocity,
and the solution of a solitary wave is expressed below the interface
displacement *η*(*x*,*t*):

in which the *η*_{0} is the amplitude. The nonlinear velocity
*V* (also called phase velocity) and the characteristic length of soliton
*L* are given as

The dispersion parameter *β* is almost larger than zero for the
internal solitary waves in the ocean, but the sign for the nonlinear
parameter *α* is relevant to the wave formation. When *α*>0, the interface displacement will show a waveform of
depression soliton. If negative, the isopycnal elevation will appear.
Therefore, the reverse situation for an internal solitary wave is determined
by the sign change of the nonlinear parameter. The KdV model is suitable
with weakly nonlinear and dispersive waves, which is capable of being used to
validate the small-amplitude ISW results in the laboratory. Nevertheless,
when nonlinearity enhancement happens because of shallower topography
or stronger stratification, the modified KdV (m-KdV) model (Michallet and
Barthélemy, 1998; Grimshaw et al., 2004) can describe relatively stronger
nonlinear solitons with the addition for the cubic nonlinearity term as

It is worth noting that the m-KdV equation takes the higher-order
nonlinear term into account and can degenerate into the KdV equation when
the cubic nonlinear parameter *α*_{1}=0. Here, the solution is
given with the interface displacement *η*(*x*,*t*):

where

More generally, when considering the continuously stratified fluid, the
linear velocity *c* refers to the long-wave velocity of each mode for the
Sturm–Liouville problem given as follows (Apel et al., 2007):

where *H* is the water depth, *N* is the buoyancy frequency, and *W* is the
nondimensional modal function when the nonlinear and dispersion parameters
(*α* and *β*, respectively) are obtained as

In addition, if still considering the background current $\stackrel{\mathrm{\u203e}}{U}\left(z\right)$, the Taylor–Goldstein equation (Miles, 1961; Liu, 2010) can
describe the vertical modal function *W* when the nonlinear and dispersion
parameters are obtained under the Boussinesq approximation expressed as
(Grimshaw et al., 2002)

where *c* is the n-mode linear speed, $\widehat{\mathit{\phi}}\left(z\right)$ is
the stream function, $\stackrel{\mathrm{\u203e}}{U}{}^{\prime \prime}$ represents the second derivative of
background currents, and *k* is the horizontal wavenumber.

This study has been conducted using the EU Copernicus Marine Service Information global ocean. Specifically, the GLORYS12V1 product in CMEMS (EU Copernicus Marine Service Information, 2022) eddy-resolving reanalysis is extracted (GLOBAL_MULTIYEAR_PHY_001_030: https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/services, last access: 25 December 2022). The current version of the nonhydrostatic ocean model (ORCTM v1) and the experiments with the internal solitary wave simulations in this paper are available through https://doi.org/10.5281/zenodo.6683597 (Huang, 2022), as are the experiment configurations, preprocessing, and post-processing. The PETSc library (The Portable, Extensible Toolkit for Scientific Computation library, https://petsc.org/release/install/download/#recommended-obtain-release-version-with-git, last access: 28 December 2022, Balay et al., 2020) needs to be installed before building the model. Nevertheless, we also provide the PETSc library of the version in use and the ORCTM quick manual for users at the above link.

HH, PS, and XC developed the nonhydrostatic dynamic framework in ORCTM and devised the internal solitary wave validation experiments. SQ and JG developed the open boundary module. HH, SQ, and XC analyzed the model results and interpreted the concepts, and all authors contributed to the writing of the paper.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank both the National Supercomputing Center in Jinan and the Marine Big Data Center of the Institute for Advanced Ocean Study at the Ocean University of China for the provision of computing resources. Also, we thank Zhigang Lai of Sun Yat-Sen University and Helmuth Haak of the Max Planck Institute for Meteorology for providing constructive comments and thoughts in model development. Finally, we would like to thank the reviewers for their careful reading of the paper and for providing professional opinions and pertinent evaluations to improve the paper.

This research has been supported by the National Key Research and Development Program of China (Super-resolution assimilation and fusion model for ocean data, SAFMOD (grant no. 2021YFF0704002)) and the National Natural Science Foundation of China (NSFC; Fundamental research on stereo-sensor detection system and seismic ocean imaging of submarine, grant no. 91958206; Study on the mechanism of generation and evolution of internal tide in the Luzon Strait influenced by background current and horizontal inhomogeneous stratification, grant no. 42276011).

This paper was edited by Xiaomeng Huang and reviewed by Andrey Pnyushkov and one anonymous referee.

Ai, C., and Ding, W.: A 3D unstructured non-hydrostatic ocean model for internal waves, Ocean Dyn., 66, 1253–1270, https://doi.org/10.1007/s10236-016-0980-9, 2016.

Ai, C., Ma, Y., Yuan, C., and Dong, G.: Non-hydrostatic model for internal wave generations and propagations using immersed boundary method, Ocean. Eng., 225, 108801, https://doi.org/10.1016/j.oceaneng.2021.108801, 2021.

Apel, J. R., Ostrovsky, L. A., Stepanyants, Y. A., and Lynch, J. F.: Internal solitons in the ocean and their effect on underwater sound, J. Acoust. Soc. Am., 121, 695–722, https://doi.org/10.1121/1.2395914, 2007.

Arakawa, A., and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, Methods. Comput. Phys., 17, 173–265, https://doi.org/10.1016/B978-0-12-460817-7.50009-4, 1977.

Arbic, B. K. and Scott, R. B.: On quadratic bottom drag, geostrophic turbulence, and oceanic mesoscale eddies, J. Phys. Oceanogr., 38, 84–103, https://doi.org/10.1175/2007JPO3653.1, 2008.

Armfield, S. and Street, R.: An analysis and comparison of the time accuracy of fractional-step methods for the Navier–Stokes equations on staggered grids, Int. J. Numer. Meth. Fl., 38, 255–282, https://doi.org/10.1002/fld.217, 2002.

Baines, P. G.: On internal tide generation models, Deep Sea. Res., 29, 307–338, https://doi.org/10.1016/0198-0149(82)90098-X, 1982.

Balay, S., Abhyankar, S., Adams, Mark F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W., Karpeyev, D., Kaushik, D., Knepley, M., May, D., McInnes, L. Curfman, Mills, R., Munson, T., Rupp, K., Sanan, P., Smith, B., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual, Argonne National Laboratory, Tech. Rep. ANL-95/11-Revision 3.13, https://doi.org/10.2172/1614847, 2020.

Benjamin, T. B.: Internal waves of finite amplitude and permanent form, J. Fluid. Mech., 25, 241–270, https://doi.org/10.1017/S0022112066001630, 1966.

Berntsen, J., Xing, J., and Alendal, G.: Assessment of non-hydrostatic ocean models using laboratory scale problems, Cont. Shelf. Res., 26, 1433–1447, https://doi.org/10.1016/j.csr.2006.02.014, 2006.

Bourgault, D. and Kelley, D. E.: A laterally averaged nonhydrostatic ocean model, J. Atmos. Ocean. Tech., 21, 1910–1924, https://doi.org/10.1175/JTECH-1674.1, 2004.

Buijsman, M. C., Kanarska, Y., and McWilliams, J. C.: On the generation and evolution of nonlinear internal waves in the South China Sea, J. Geophys. Res.-Oceans, 115, C02012, https://doi.org/10.1029/2009JC005275, 2010a.

Buijsman, M. C., McWilliams, J. C., and Jackson, C. R.: East-west asymmetry in nonlinear internal waves from Luzon Strait, J. Geophys. Res.-Oceans, 115, C10057, https://doi.org/10.1029/2009JC006004, 2010b.

Chen, C., Liu, H., and Beardsley, R. C.: An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: application to coastal ocean and estuaries, J. Atmos. Ocean. Tech., 20, 159–186, https://doi.org/10.1175/1520-0426(2003)020<0159:AUGFVT>2.0.CO;2, 2003.

Chen, X., Jungclaus, J., Thomas, M., Maier-Reimer, E., Haak, H., and Suendermann, J.: An oceanic general circulation and tide model in orthogonal curvilinear coordinates, Amer. Geophys. Union., Fall Meeting 2005, San Francisco, CA, December 2005, Abstract OS41B-0600, https://ui.adsabs.harvard.edu/abs/2005AGUFMOS41B0600C/abstract (last access: 25 December 2022), 2005.

Chen, Z., Nie, Y., Xie, J., Xu, J., He, Y., and Cai, S.: Generation of internal solitary waves over a large sill: From Knight Inlet to Luzon Strait, J. Geophys. Res.-Oceans, 122, 1555–1573. https://doi.org/10.1002/2016JC012206, 2017.

Chorin, A. J.: Numerical solution of the Navier–Stokes equations, Math. Comput., 22, 745–762, https://doi.org/10.2307/2004575, 1968.

Cushman-Roisin, B.: Kelvin–Helmholtz instability as a boundary-value problem, Environ. Fluid. Mech., 5, 507–525, https://doi.org/10.1007/s10652-005-2234-0, 2005.

Duda, T. F., Morozov, A. K., Howe, B. M., Brown, M. G., Speer, K., Lazarevich, P., Worcester, P. F., and Cornuelle, B. D.: Evaluation of a long-range joint acoustic navigation/thermometry system, Oceans 2006 IEEE, 1–6, https://doi.org/10.1109/OCEANS.2006.306999, 2006.

EU Copernicus Marine Service Information: Global Ocean Physics Reanalysis: GLOBAL_MULTIYEAR_PHY_001_030, Copernicus.eu [data set], https://doi.org/10.48670/moi-00021, 2022.

Fofonoff, N. P. and Millard Jr., R. C.: Algorithms for computation of fundamental properties of seawater, Paris, France, UNESCO, 53 pp., https://doi.org/10.25607/OBP-1450, 1983.

Fringer, O. B., Gerritsen, M., and Street, R. L.: An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator, Ocean. Model., 14, 139–173, https://doi.org/10.1016/j.ocemod.2006.03.006, 2006.

Garrett, C. and Kunze, E.: Internal tide generation in the deep ocean, Annu. Rev. Fluid. Mech., 39, 57–87, https://doi.org/10.1146/annurev.fluid.39.050905.110227, 2007.

Gerkema, T., and Zimmerman, J. T. F.: Generation of nonlinear internal tides and solitary waves, J. Phys. Oceanogr., 25, 1081–1094, https://doi.org/10.1175/1520-0485(1995)025<1081:GONITA>2.0.CO;2, 1995.

Gilbert, D., and Garrett, C.: Implications for ocean mixing of internal wave scattering off irregular topography, J. Phys. Oceanogr., 19, 1716–1729, https://doi.org/10.1175/1520-0485(1989)019<1716:IFOMOI>2.0.CO;2, 1989.

Grimshaw, R., Pelinovsky, E., and Poloukhina, O.: Higher-order Korteweg-de Vries models for internal solitary waves in a stratified shear flow with a free surface, Nonlin. Processes Geophys., 9, 221–235, https://doi.org/10.5194/npg-9-221-2002, 2002.

Grimshaw, R., Pelinovsky, E., Talipova, T., and Kurkin, A.: Simulation of the transformation of internal solitary waves on oceanic shelves, J. Phys. Oceanogr., 34, 2774–2791, https://doi.org/10.1175/JPO2652.1, 2004.

Grue, J., Jensen, A., Rusås, P. O., and Sveen, J. K.: Breaking and broadening of internal solitary waves, J. Fluid. Mech., 413, 181–217, https://doi.org/10.1017/S0022112000008648, 2000.

Härtel, C., Meiburg, E., and Necker, F.: Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries, J. Fluid. Mech., 418, 189–212, https://doi.org/10.1017/s0022112000001221, 2000.

Huang, H.: HuangOCEAN02/ORCTM: ORCTM v1.1.1 (ORCTMv1.1.1), Zenodo [code], https://doi.org/10.5281/zenodo.6683597, 2022.

Helfrich, K. R. and Melville, W. K.: On long nonlinear internal waves over slope-shelf topography, J. Fluid. Mech., 167, 285–308, https://doi.org/10.1017/S0022112086002823, 1986.

Huang, X., Chen, Z., Zhao, W., Zhang, Z., Zhou, C., Yang, Q., and Tian, J.: An extreme internal solitary wave event observed in the northern South China Sea, Sci. Rep.-UK, 6, 1–10, https://doi.org/10.1038/srep30041, 2016.

Kanarska, Y., Shchepetkin, A., and McWilliams, J. C.: Algorithm for non-hydrostatic dynamics in the regional oceanic modeling system, Ocean. Model., 18, 143–174, https://doi.org/10.1016/j.ocemod.2007.04.001, 2007.

Ko, D. S., Martin, P. J., Rowley, C. D., and Preller, R. H.: A real-time coastal ocean prediction experiment for MREA04, J. Marine. Syst., 69, 17–28, https://doi.org/10.1016/j.jmarsys.2007.02.022, 2008.

Lai, Z., Chen, C., Cowles, G. W., and Beardsley, R. C.: A nonhydrostatic version of FVCOM: 1. Validation experiments, J. Geophys. Res.-Oceans, 115, C11010, https://doi.org/10.1029/2009JC005525, 2010.

Lawrence, G. A., Browand, F. K., and Redekopp, L. G.: The stability of a sheared density interface, Phys. Fluids. A-Fluid., 3, 2360–2370, https://doi.org/10.1063/1.858175, 1991.

Legg, S., and Adcroft, A.: Internal wave breaking at concave and convex continental slopes, J. Phys. Oceanogr., 33, 2224–2246, https://doi.org/10.1175/1520-0485(2003)033<2224:IWBACA>2.0.CO;2, 2003.

Legg, S. and Klymak, J.: Internal hydraulic jumps and overturning generated by tidal flow over a tall steep ridge, J. Phys. Oceanogr., 38, 1949–1964, https://doi.org/10.1175/2008JPO3777.1, 2008.

Li, Q.: Numerical assessment of factors affecting nonlinear internal waves in the South China Sea, Prog. Oceanogr., 121, 24–43, https://doi.org/10.1016/j.pocean.2013.03.006, 2014.

Li, J., Zhang, Q., and Chen, T.: ISWFoam: a numerical model for internal solitary wave simulation in continuously stratified fluids, Geosci. Model Dev., 15, 105–127, https://doi.org/10.5194/gmd-15-105-2022, 2022.

Liu, Z.: Instability of baroclinic tidal flow in a stratified fjord, J. Phys. Oceanogr., 40, 139–154, https://doi.org/10.1175/2009JPO4154.1, 2010.

Liu, Z., Lin, L., Xie, L., and Gao, H.: Partially implicit finite difference scheme for calculating dynamic pressure in a terrain-following coordinate non-hydrostatic ocean model, Ocean. Model., 106, 44–57, https://doi.org/10.1016/j.ocemod.2016.09.004, 2016.

Ma, Q., Yuan, C., Lin, X., and Chen, X.: The investigation of internal solitary waves over a continental shelf-slope, J. Oceanol. Limnol., 38, 695–706, https://doi.org/10.1007/s00343-019-9123-8, 2020.

Marshall, J., Hill, C., Perelman, L., and Adcroft, A.: Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling, J. Geophys. Res.-Oceans, 102, 5733–5752, https://doi.org/10.1029/96JC02776, 1997a.

Marshall, J., Adcroft, A., Hill, C., Perelman L., and Heisey, C.: A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers, J. Geophys. Res.-Oceans, 102, 5753–5766, https://doi.org/10.1029/96JC02775, 1997b.

Marshall, J., Jones, H., and Hill, C.: Efficient ocean modeling using non-hydrostatic algorithms, J. Marine. Syst., 18, 115–134, https://doi.org/10.1016/S0924-7963(98)00008-6, 1998.

Marsland, S. J., Haak, H., Jungclaus, J. H., Latif, M. and Röske, F.: The Max-Planck-Institute global ocean/sea ice model with orthogonal curvilinear coordinates, Ocean. Model., 5, 91–127, https://doi.org/10.1016/S1463-5003(02)00015-X, 2003.

Michallet, H. and Barthélemy, E.: Experimental study of interfacial solitary waves, J. Fluid. Mech., 366, 159–177, https://doi.org/10.1017/S002211209800127X, 1998.

Michallet, H. and Ivey, G. N.: Experiments on mixing due to internal solitary waves breaking on uniform slopes, J. Geophys. Res.-Oceans, 104, 13467–13477, https://doi.org/10.1029/1999JC900037, 1999.

Miles, J. W.: On the stability of heterogeneous shear flows, J. Fluid. Mech., 10, 496–508, https://doi.org/10.1017/S0022112061000305, 1961.

Mtfller, P.: On the diffusion of momentum and mass by internal gravity waves, J. Fluid. Mech., 77, 789–823, https://doi.org/10.1017/S0022112076002899, 1976.

Ono, H.: Algebraic solitary waves in stratified fluids, J. Phys. Soc. Jpn., 39, 1082–1091, https://doi.org/10.1143/JPSJ.39.1082, 1975.

Osborne, A. R., Burch, T. L., and Scarlet, R. I.: The influence of internal waves on deep-water drilling, J. Pet. Technol., 30, 1497–1504, https://doi.org/10.2118/6913-PA, 1978.

Pacanowski, R. C., and Philander, S. G. H.: Parameterization of vertical mixing in numerical models of tropical oceans, J. Phys Oceanogr., 11, 1443-1451, https://doi.org/10.1175/1520-0485(1981)011<1443:POVMIN>2.0.CO;2, 1981.

Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T.: Numerical Recipes in C. The Art of Scientific Computing, Cambridge University Press, ISBN: 0-521-35465-X, 1988.

Ramp, S. R., Tang, T. Y., Duda, T. F., Lynch, J. F., Liu, A. K., Chiu, C. S., Bahr, F. L., Kim, H. R., and Yang, Y. J.: Internal solitons in the northeastern South China Sea. Part I: Sources and deep water propagation, IEEE. J. Oceanic. Eng., 29, 1157–1181, https://doi.org/10.1109/JOE.2004.840839, 2004.

Ramp, S. R., Park, J. -H., Yang, Y. J., Bahr, F. L., and Jeon, C.: Latitudinal Structure of Solitons in the South China Sea, J. Phys. Oceanogr., 49, 1747–1767, https://doi.org/10.1175/JPO-D-18-0071.1, 2019.

Saad Y.: A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Comput., 14, 461–469, https://doi.org/10.1137/0914028, 1993.

Shaw, P. T., Ko, D. S., and Chao, S. Y.: Internal solitary waves induced by flow over a ridge: With applications to the northern South China Sea, J. Geophys. Res.-Oceans, 114, C02019, https://doi.org/10.1029/2008JC005007, 2009.

Shchepetkin, A. F., and McWilliams, J. C.: The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model, Ocean. Model., 9, 347–404, https://doi.org/10.1016/j.ocemod.2004.08.002, 2005.

Smith, B., Bjørstad, P., and Gropp, W.: Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, https://doi.org/10.1016/S0898-1221(97)90035-3, 1996.

Stansby, P. K. and Zhou, J. G.: Shallow‐water flow solver with non‐hydrostatic pressure: 2D vertical plane problems, Int. J. Numer. Meth. Fluids, 28, 541–563, https://doi.org/10.1002/(SICI)1097-0363(19980915)28:3<541::AID-FLD738>3.0.CO;2-0, 1998.

Vlasenko, V. and Hutter, K.: Numerical experiments on the breaking of solitary internal waves over a slope–shelf topography, J. Phys. Oceannogr., 32, 1779–1793, https://doi.org/10.1175/1520-0485(2002)032<1779:NEOTBO>2.0.CO;2, 2002.

Vlasenko, V., Stashchuk, N., and Hutter, K.: Baroclinic tides: theoretical modeling and observational evidence, Cambridge University Press, ISBN: 978-0-521-84395-9, 2005.

Vlasenko, V., Stashchuk, N., Guo, C., and Chen, X.: Multimodal structure of baroclinic tides in the South China Sea, Nonlin. Processes Geophys., 17, 529–543, https://doi.org/10.5194/npg-17-529-2010, 2010.

Wang, Y. H., Dai, C. F., and Chen, Y. Y.: Physical and ecological processes of internal waves on an isolated reef ecosystem in the South China Sea, Geophys. Res. Lett., 34, 312–321, https://doi.org/10.1029/2007GL030658, 2007.

Wessels, F. and Hutter, K.: Interaction of internal waves with a topographic sill in a two-layered fluid, J. Phys. Oceanogr., 26, 5–20, https://doi.org/10.1175/1520-0485(1996)026<0005:IOIWWA>2.0.CO;2, 1996.

Wolff, J. -O., Maier-Reimer, E., and Legutke, S.: The Hamburg ocean primitive equation model, Tech. Rep. No. 13, German Climate Computer Center (DKRZ), Hamburg, Germany, ISSN: 0940-9327, https://hdl.handle.net/21.11116/0000-0009-1487-8 (last access: 25 December 2022), 1997.

Zhang, Z., Fringer, O. B., and Ramp, S. R.: Three-dimensional, nonhydrostatic numerical simulation of nonlinear internal wave generation and propagation in the South China Sea, J. Geophys. Res.-Oceans, 116, C05022, https://doi.org/10.1029/2010JC006424, 2011.

Zhao, Z. and Alford, M. H.: Source and propagation of internal solitary waves in the northeastern South China Sea, J. Geophys. Res.-Oceans, 111, C11012, https://doi.org/10.1029/2006JC003644, 2006.

Zeng, Z., Chen, X., Yuan, C., Tang, S., and Chi, L.: A numerical study of generation and propagation of type-a and type-b internal solitary waves in the northern South China Sea, Acta Oceanol. Sin., 38, 20–30, https://doi.org/10.1007/s13131-019-1495-2, 2019.

Zheng, Q., Susanto, R. D., Ho, C. R., Song, Y. T., and Xu, Q.: Statistical and dynamical analyses of generation mechanisms of solitary internal waves in the northern South China Sea, J. Geophys. Res.-Oceans, 112, C0302, https://doi.org/10.1029/2006JC003551, 2007.

- Abstract
- Introduction
- Model development
- Model applications and assessments
- Discussion and conclusion
- Appendix A: Discretization algorithms of the Poisson equation
- Appendix B: The Korteweg–de Vries (KdV) model in shallow water
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Model development
- Model applications and assessments
- Discussion and conclusion
- Appendix A: Discretization algorithms of the Poisson equation
- Appendix B: The Korteweg–de Vries (KdV) model in shallow water
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References