ISWFoam: a numerical model for internal solitary wave simulation in continuously stratified fluids

A numerical model, ISWFoam, for simulating internal solitary waves (ISWs) in continuously stratified, incompressible, viscous fluids is developed based on a fully three-dimensional (3D) Navier–Stokes equation using the open-source code OpenFOAM^®. This model combines the density transport equation with the Reynolds-averaged Navier–Stokes equation with the Coriolis force, and the model discrete equation adopts the finite-volume method. The k–ω SST turbulence model has also been modified according to the variable density field. ISWFoam provides two initial wave generation methods to generate an ISW in continuously stratified fluids, including solving the weakly nonlinear models of the extended Korteweg–de Vries (eKdV) equation and the fully nonlinear models of the Dubreil–Jacotin–Long (DJL) equation. Grid independence tests for ISWFoam are performed, and considering the accuracy and computing efficiency, the appropriate grid size of the ISW simulation is recommended to be 1$/$150th of the characteristic length and 1$/$25th of the ISW amplitude. Model verifications are conducted through comparisons between the simulated and experimental data for ISW propagation examples over a flat bottom section, including laboratory scale and actual ocean scale, a submerged triangular ridge, a Gaussian ridge, and slope. The laboratory test results, including the ISW profile, wave breaking location, ISW arrival time, and the spatial and temporal changes in the mixture region, are well reproduced by ISWFoam. The ISWFoam model with unstructured grids and local mesh refinement can effectively simulate the evolution of ISWs, the ISW breaking phenomenon, waveform inversion of ISWs, and the interaction between ISWs and complex topography.

1 Introduction

Internal solitary waves (ISWs) are commonly observed in oceans, particularly in continental shelf regions, due to strong tidal current flows over large topographic features (Huthnance, 1981), such as in the northern South China Sea (Alford et al., 2010, 2015; Cai et al., 2012). ISWs play an important role in both conveying nutrients from the deep ocean to shallower layers and promoting biological growth (Sandstrom et al., 1984). Additionally, ISWs are a potential threat to the ocean structures of resource exploration, exploitation and submarine navigation vehicles (Alford et al., 2010; Osborne and Burch, 1980). A considerable number of studies, which include field measurements, remote sensing, experiments, theoretical analysis and numerical simulations, have been carried out due to the significance of ISWs (Vlasenko et al., 2005; Apel et al., 2006; Alford et al., 2011; Guo et al., 2014).

For numerically simulated ISWs, many models have been adopted, including the Euler equation, the inviscid or viscid incompressible Boussinesq model, the hydrostatic model, the non-hydrostatic model, and the VOF-based two-phase flow model. Among these models, the representative hydrostatic models include the Naval Research Laboratory Ocean Nowcast/Forecast System (ONFS) (Ko et al., 2008), the Regional Hallberg Isopycnal Tide Mode (RHIMT) (Hallberg and Rhines, 1996; Hallberg, 1997) and the Ostrovsky–Hunter model. The representative non-hydrostatic models include the Bergen Ocean Model (BOM), the nonhydrostatic Regional Ocean Modeling System model (ROMS), the Stanford Unstructured Nonhydrostatic Terrain-following Adaptive Navier–Stokes Simulator (SUNTANS) and the Massachusetts Institute of Technology general circulation model (MITgcm). For example, Zhang et al. (2012) established a variable-water-depth internal-wave numerical model in a continuously stratified fluid system based on the Euler equation. Xu and Stastna (2020) used the viscid incompressible Boussinesq model to study cross-boundary-layer transport (Boegman and Stastna, 2019) by the fissioning process of shoaling ISWs. Lamb (1994) established a non-hydrostatic model, using a second-order projection method developed by Bell and Marcus (1992), which is used for internal-wave research including boundary-layer instability (Aghsaee et al., 2012), reflection (Lamb and Nguyen, 2009 ), and the interaction of the tides with the topography (Lamb, 2007; Aghsaee et al., 2010). Diamessis (2005) developed a spectral multidomain penalty method model and correctly reproduced the characteristic vorticity and internal-wave structure. Subich et al. (2013) developed a spectral collocation method for the solution of the Navier–Stokes equations under the Boussinesq approximation and simulated the internal wave in continuously stratified fluid. Smedstad et al. (2003) employed the ONFS model to establish a global ocean real-time forecasting system with an operational eddy resolution of $\mathrm{1}/\mathrm{16}$^∘, which effectively tracks ocean eddies, ocean currents and ocean fronts. Simmons et al. (2004) employed the RHIMT model to carry out a global numerical simulation of tidal currents and analysed the whole process of the conversion rate of barotropic waves into baroclinic waves. Thiem (2011) used the Bergen Ocean Model to explore the bottom boundary-layer flow caused by waves beneath a propagating ISW in a two-fluid system. Li and Farmer (2011) employed the Ostrovsky–Hunter model to study the nonlinear evolution of a monochromatic internal wave. Buijsman et al. (2010) employed a ROMS model to study the asymmetry in solitons to the east and west of the Luzon Strait. Zhang et al. (2011) used the nonhydrostatic SUNTANS model (Fringer et al., 2006) to study the dynamics of A-wave and B-wave formation. Rayson et al. (2018) used the modified SUNTANS model to study the internal waves around Scott Reef and provided the generation process of internal lee waves. Vlasenko et al. (2010) employed the MITgcm model to investigate the baroclinic tidal energy conversion in the area west of the Luzon Strait.

In summary, for continuously stratified fluids in complex ocean environments, numerical simulation has become a leading method for ISW investigations. However, there are presently few versatile numerical models with shared code that can accurately simulate the ISW flow around complex topography and submarine navigation vehicles in continuously stratified fluids. Therefore, the main objective of this paper is to develop a solver, referred to as ISWFoam with a modified k–ω SST model that considers the variable density field, which simulates the ISW in continuously stratified, incompressible and viscous fluids using the finite-volume method with unstructured grids based on a fully three-dimensional (3D) Navier–Stokes equation using the OpenFOAM^® library.

Notably, the open-source field operation and manipulation code OpenFOAM, as an object-oriented C$++$ open-source library that can be used to build a variety of solvers for computational fluid problems based on the finite-volume method, is becoming increasingly popular in the computational fluid research community. At present, the official version of OpenFOAM does not have a solver or boundary conditions for solving the ISW in continuously stratified fluids. Although some researchers simulate ISWs by modifying the OpenFOAM code, most of these studies are based on a two-fluid system without considering continuous stratification in density, such as Meng and Zhang (2016) and Li et al. (2017). Though recent work by Ding et al. (2020) and Li et al. (2021) considered continuous stratification in density, their wave generation theories do not consider continuous stratification in density. To extensively use the numerical model of ISWs as a tool in the future, we will develop ISWFoam to simulate the ISW in continuously stratified, incompressible and viscous fluids based on the OpenFOAM library. The turbulence model will consider the variable density field. In addition, ISWFoam will provide two initial methods to generate an ISW in continuously stratified fluids, including solving the weakly nonlinear models of the extended Korteweg–de Vries (eKdV) equation and the fully nonlinear models of the Dubreil–Jacotin–Long (DJL) equation. This approach renders the numerical model suitable for the simulation of ISW flows in complex geometries and topographies. It is worth noting that ISWFoam does not consider the generation process of ISWs but focuses on the propagation and evolution of ISWs that have already been generated, and the interaction between ISWs and complex structures and topography on field scales.

The outline of the paper is described as follows. First, in Sect. 2, the governing equations for a continuously stratified fluid are presented, and discrete forms of these equations are derived. Then, grid independence tests of the developed ISWFoam model are described in Sect. 3. Subsequently, in Sect. 4, a series of test cases are presented to verify the model. Simulation examples at the field scale are shown in Sect. 5. Finally, the conclusions are drawn in Sect. 6.

2 ISWFoam: A three-dimensional numerical solver for ISWs in a continuously stratified fluid

2.1 Governing equations

We present an ISW numerical model by solving the motion of a three-dimensional, viscous, incompressible fluid with the Boussinesq approximation and rigid lid hypothesis. The governing equations of the model are

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{U}=\mathrm{0},\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{U}}{\partial t}}& +\left(\mathit{U}\cdot \mathrm{\nabla }\right)\mathit{U}-\mathrm{\nabla }\cdot \left({\mathit{\nu }}_{\mathrm{Eff}}\mathrm{\nabla }\mathit{U}\right)=Q\\ Q& ={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\left(-\mathrm{\nabla }{p}_{\mathit{\_}\mathrm{rgh}}-\mathit{g}\cdot \mathit{X}\mathrm{\nabla }\mathit{\rho }-\mathrm{\Omega }{\mathit{e}}_{\mathrm{3}}\times \mathit{U}\right),\end{array}\text{(3)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\left(\mathit{U}\cdot \mathrm{\nabla }\right)\mathit{\rho }=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }\mathit{\rho }\right),\end{array}$$

where U = (u_i, u_j, u_k) is the velocity vector, t is time, ∇ is the gradient operator, Q is the source term, ρ₀ is the reference density, ρ is the density field, p_{_rgh} is a modified pressure field, g is the gravitational acceleration vector and X is the position vector. v_Eff is the effective kinematic viscosity defined as $v_{\mathrm{Eff}}={\mathit{\mu }}_{\mathrm{Eff}}/{\mathit{\rho }}_{\mathrm{0}}$, where μ_Eff is the effective dynamic viscosity including the molecular viscosity (μ_l) and turbulent viscosity (μ_t). k is the diffusion coefficient, and its value is the same as the effective dynamic viscosity(μ_Eff). Ω is the Coriolis parameter, which is twice the speed of rotation around the vertical unit vector e₃= (0, 0, 1). ISWFoam uses a modified pressure p_{_rgh} instead of a total pressure p, and their relationship is given by

$$\begin{array}{}\text{(4)}& p_{\mathit{\_}\mathrm{rgh}}=p-\mathit{\rho }\mathit{g}\cdot \mathit{X},\mathrm{\nabla }{p}_{\mathit{\_}\mathrm{rgh}}=\mathrm{\nabla }p-\mathit{\rho }\mathit{g}-\mathit{g}\cdot \mathit{X}\mathrm{\nabla }\mathit{\rho },\end{array}$$

The upper boundary (z=H, with H the depth of computation domain) is treated as a rigid lid, and the kinematic boundary conditions for this boundary are given by

$$\begin{array}{}\text{(5)}& u_k(x,y,H,t)=\mathrm{0}.\end{array}$$

To close the above equations, the turbulence model needs to be employed. The two-equation k–ε model is widely used as an effective turbulence model, but it cannot capture the proper behaviour of turbulent boundary layers up to separation due to adverse pressure gradients (Wilcox, 1993). For the above boundary-layer separation problem, Bardina et al. (1997) and Menter et al. (2003) suggested the use of the k–ω shear stress transport (SST) model to obtain substantially more accurate results. Therefore, the turbulence model used in this paper is the k–ω SST model. Note that in OpenFOAM, the incompressible version for turbulence models does not consider the variable density field, and instead, it treats the density as a constant, such as in the k–ω SST model:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial k}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{U}k\right)=\mathrm{\nabla }\cdot \left[\left({\mathit{\nu }}_{\mathrm{Eff}}+{\mathit{\sigma }}_{\mathrm{k}}{\mathit{\nu }}_t\right)\mathrm{\nabla }k\right]+P_k^{\ast }-{\mathit{\beta }}^{\ast }\mathit{\omega }k,\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{\omega }}{\partial t}}& +\mathrm{\nabla }\cdot \left(\mathit{U}\mathit{\omega }\right)=\mathrm{\nabla }\cdot \left[\left({\mathit{\nu }}_{\mathrm{Eff}}+{\mathit{\sigma }}_{\mathit{\omega }}{\mathit{\nu }}_t\right)\mathrm{\nabla }\mathit{\omega }\right],\\ & +C_{\mathit{\gamma }}{\displaystyle \frac{\mathit{\omega }}{k}}{P}_k-C_{\mathit{\beta }}{\mathit{\omega }}^{\mathrm{2}}+\mathrm{2}\left(\mathrm{1}-F_{\mathrm{1}}\right){\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\omega }\mathrm{2}}}{\mathit{\omega }}}\mathrm{\nabla }k\cdot \mathrm{\nabla }\mathit{\omega },\end{array}\text{(8)}& {\displaystyle }{P}_k^{\ast }=\min (P_k,c_{\mathrm{1}}{C}_{\mathit{\mu }}k\mathit{\omega }),\text{(9)}& {\displaystyle }{\mathit{\nu }}_t={\displaystyle \frac{a_{\mathrm{1}}k}{max\left(a_{\mathrm{1}}\mathit{\omega },\sqrt{\mathrm{2}}{S}_t{F}_{\mathrm{2}}\right)}},\end{array}$$

where k is the turbulent kinetic energy, ω is the specific dissipation rate; P_k is the production term of k; $P_k={\mathit{\tau }}^R:\mathrm{\nabla }\mathit{U}$; $P_k^{\ast }$ is related to the production term of turbulence kinetic energy P_k in the k equation; v_t is the turbulent kinematic viscosity; S_t is the mean rate of the flow strain; S_t=0.5(∇U+∇U^T); the model constants are assigned the values ${\mathit{\beta }}^{\ast }=\mathrm{0.09}$, a₁=0.31, c₁=10 and C_μ=0.09; F₁ and F₂ are blending functions; and the value of σ_k, σ_ω, C_γ and C_β are blended using the equation $\mathit{Ô}=F_{\mathrm{1}}{\mathit{Ô}}_{\mathrm{1}}+(\mathrm{1}-F_{\mathrm{1}}){\mathit{Ô}}_{\mathrm{2}}$ in which Ô₁ and Ô₂ are given in Table 1.

Table 1 Default values for Ô₁ and Ô₂.

Considering the variable density field during the solution process, it is necessary to consider the change in the density field in the turbulence model. Therefore, we modify the turbulence model to consider the change in density, and finally a modified k–ω SST model that considers the change in density is used to close the equation:

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{\rho }k}{\partial t}}& +\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{U}k\right)=\mathrm{\nabla }\cdot \left[\mathit{\rho }\left({\mathit{\nu }}_{\mathrm{Eff}}+{\mathit{\sigma }}_k{\mathit{\nu }}_t\right)\mathrm{\nabla }k\right]\\ & +\mathit{\rho }{P}_k^{\ast }-\mathit{\rho }{\mathit{\beta }}^{\ast }\mathit{\omega }k\end{array}\text{(11)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{\rho }\mathit{\omega }}{\partial t}}& +\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{U}\mathit{\omega }\right)=\mathrm{\nabla }\cdot \left[\mathit{\rho }\left({\mathit{\nu }}_{\mathrm{Eff}}+{\mathit{\sigma }}_{\mathit{\omega }}{\mathit{\nu }}_t\right)\mathrm{\nabla }\mathit{\omega }\right]\\ & +C_{\mathit{\gamma }}{\displaystyle \frac{\mathit{\omega }}{k}}{P}_k-C_{\mathit{\beta }}\mathit{\rho }{\mathit{\omega }}^{\mathrm{2}}+\mathrm{2}\left(\mathrm{1}-F_{\mathrm{1}}\right)\mathit{\rho }{\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\omega }\mathrm{2}}}{\mathit{\omega }}}\mathrm{\nabla }k\cdot \mathrm{\nabla }\mathit{\omega }.\end{array}\end{array}$$

2.2 Numerical discretisation

The governing equations are numerically discretised using the finite-volume method based on the C$++$ open-source library of OpenFOAM. The finite-volume method requires that Eqs. (2) and (3) are satisfied over the control volume V_P around point P in integral form:

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}\underset{V_P}{\int }\underset{\mathrm{\Delta }t}{\int }& \left[{\displaystyle \frac{\partial \mathit{U}}{\partial t}}+\left(\mathit{U}\cdot \mathrm{\nabla }\right)\mathit{U}-\mathrm{\nabla }\cdot \left({\mathit{\nu }}_{\mathrm{Eff}}\mathrm{\nabla }\mathit{U}\right)\right]\mathrm{d}V\mathrm{d}t\\ & =\underset{V_P}{\int }\underset{\mathrm{\Delta }t}{\int }Q\mathrm{d}V\mathrm{d}t,\end{array}\text{(13)}& {\displaystyle }\underset{V_P}{\int }\underset{\mathrm{\Delta }t}{\int }\left[{\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\left(\mathit{U}\cdot \mathrm{\nabla }\right)\mathit{\rho }-\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }\mathit{\rho }\right)\right]\mathrm{d}V\mathrm{d}t=\mathrm{0},\end{array}$$

where Δt is the time step.

The momentum equation in ISWFoam is solved by constructing a predicted velocity field and then using the pressure implicit with splitting of operators (PISO) algorithm (Issa, 1986) to modify it. n is defined to represent the current moment. The PISO iteration process is marked as m; when m is equal to zero, it represents the current time (tⁿ).

First, only the temporal, convection and diffusion terms appear in the discrete version of the equation momentum, and the other terms are ignored. After this operation, we obtain an explicit expression for the predicted velocity field ${\mathit{U}}_P^r$, namely,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\displaystyle \frac{{\mathit{U}}_P^r-{\mathit{U}}_P^n}{\mathrm{\Delta }t}}{V}_P& +\sum _{f\in \partial V_P}\left({\mathit{\varphi }}_f^n{\mathit{U}}_{\mathrm{f}}^r\right)\\ & -\sum _{f\in \partial V_P}{\mathit{\nu }}_{\mathrm{Eff}}\mathrm{\nabla }{\mathit{U}}_{\mathrm{f}}^r\cdot {\mathit{S}}_{\mathrm{f}}=\mathrm{0},\end{array}\end{array}$$

where P represents the centre of the grid cell, ${\mathit{\varphi }}_f^n={\mathit{U}}_{\mathrm{f}}^n\cdot {\mathit{S}}_{\mathrm{f}}$ is the volume flux at the initial time n and S_f is the face vector.

The solution process requires the velocity on the surface f. Assuming the variation in ${\mathit{U}}_{\mathrm{f}}^r$ between the centre P of the grid and the centre N of the adjacent grid, the face values are calculated using a mixture method (blended differencing) of the central scheme (central differencing) and the upwind scheme (upwind differencing) as follows (Jasak, 1996):

$$\begin{array}{}\text{(15)}& {\mathit{U}}_{\mathrm{f}}=\left(\mathrm{1}-{\mathit{\lambda }}_U\right){\left({\mathit{U}}_{\mathrm{f}}\right)}_{UD}+{\mathit{\lambda }}_U{\left({\mathit{U}}_{\mathrm{f}}\right)}_{CD},\end{array}$$

where

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\left({\mathit{U}}_{\mathrm{f}}\right)}_{UD}& =\left\{\begin{array}{l}{\mathit{U}}_P\mathrm{for}{\mathit{\varphi }}_{\mathrm{f}}\ge \mathrm{0},\\ {\mathit{U}}_N\mathrm{for}{\mathit{\varphi }}_{\mathrm{f}}<\mathrm{0},\end{array}\right.\mathrm{and}{\left({\mathit{U}}_{\mathrm{f}}\right)}_{CD}\\ & ={\displaystyle \frac{{\mathit{U}}_P+{\mathit{U}}_N}{\mathrm{2}}},\end{array}\end{array}$$

where N represents the centre of the adjacent grid cells, and ${\mathit{\varphi }}_{\mathrm{f}}={\mathit{U}}_{\mathrm{f}}\cdot {\mathit{S}}_{\mathrm{f}}$ is volume flux. The limiter λ_U can be selected from several alternatives (OpenFOAM, 2019), including linear, QUICK, van Leer, etc. In the following derivation process, the van Leer scheme was used to calculate the velocity of the face centre

$$\begin{array}{}\text{(17)}& {\mathit{U}}_{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{U}}_P+{\mathit{U}}_N\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathit{\psi }\left({\mathit{\varphi }}_f\right)(\mathrm{1}-{\mathit{\lambda }}_U)\right]\left({\mathit{U}}_P-{\mathit{U}}_N\right),\end{array}$$

where ψ(ϕ_f) is a step function defined by

$$\begin{array}{}\text{(18)}& \mathit{\psi }\left({\mathit{\varphi }}_f\right)=\left\{\begin{array}{l}\mathrm{1}\mathrm{for}{\mathit{\varphi }}_f\ge \mathrm{0},\\ -\mathrm{1}\mathrm{for}{\mathit{\varphi }}_f<\mathrm{0}.\end{array}\right.\end{array}$$

Inserting Eq. (15) into Eq. (13) yields

$$\begin{array}{}\text{(19)}& A_P{\mathit{U}}_P^r=\sum _{f\in \partial V_P}{A}_N{\mathit{U}}_P^m+{\displaystyle \frac{{\mathit{U}}_P^n}{\mathrm{\Delta }t}}=H\left({\mathit{U}}^m\right).\end{array}$$

After some manipulation, the quantities A_P and A_N are given as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{A}_P=\left\{{\displaystyle \frac{V_P}{\mathrm{\Delta }t}}\right.& +\sum _{f\in \partial V_P}{\displaystyle \frac{{\mathit{\varphi }}_f^n}{\mathrm{2}}}\left[\mathrm{1}+\mathit{\psi }\left({\mathit{\varphi }}_f\right)(\mathrm{1}-{\mathit{\lambda }}_U)\right]\\ & +\sum _{f\in \partial V_P}{\mathit{\nu }}_{\mathrm{Eff},f}{\displaystyle \frac{\left|{\mathit{S}}_{\mathrm{f}}\right|}{\left|d\right|}}{\displaystyle \frac{\mathrm{1}}{V_P}},\end{array}\end{array}$$

$$\begin{array}{}\text{(21)}& A_N=\left\{-{\displaystyle \frac{{\mathit{\varphi }}_f^n}{\mathrm{2}}}\left[\mathrm{1}-\mathit{\psi }\left({\mathit{\varphi }}_f\right)(\mathrm{1}-{\mathit{\lambda }}_U)\right]+{\mathit{\nu }}_{\mathrm{Eff},f}{\displaystyle \frac{\left|{\mathit{S}}_{\mathrm{f}}\right|}{\left|d\right|}}{\displaystyle \frac{\mathrm{1}}{V_P}}.\right.\end{array}$$

Including the effect of gravity and the Coriolis force in Eq. (17),

$$\begin{array}{}\text{(22)}& {\mathit{U}}_P^r={\displaystyle \frac{H\left({\mathit{U}}^m\right)}{A_P}}-{\displaystyle \frac{{\left(\mathit{g}\cdot \mathit{X}\mathrm{\nabla }\mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}\right)}^n}{A_P}}-{\displaystyle \frac{{\left(\mathrm{\Omega }{\mathit{e}}_{\mathrm{3}}\times \mathit{U}\right)}^n}{A_P}}.\end{array}$$

Notably, that when m is equal to zero, it represents the initial moment n, and the value of the initial moment is known. Therefore, we obtain the predicted velocity field ${\mathit{U}}_P^r$ in the first iteration. We define the surface gradient operator $\left(\mathrm{\nabla }\frac{\mathrm{1}}{f}\right)$, and the type of gradient operator acting on U is $\mathrm{\nabla }\frac{\mathrm{1}}{f}\mathit{U}=\left({\mathit{U}}_N^m-{\mathit{U}}_P^m\right)/\left|d\right|$, which represents the distance from the centre of the grid N to P. Similarly, the surface gradient operator $\left(\mathrm{\nabla }\frac{\mathrm{1}}{f}\right)$ acting on scalar γ is $\mathrm{\nabla }\frac{\mathrm{1}}{f}\mathit{\lambda }=\left({\mathit{\lambda }}_N^m-{\mathit{\lambda }}_P^m\right)/\left|d\right|$. The associated flux (${\mathit{\varphi }}_f={\mathit{U}}_{\mathrm{f}}\cdot {\mathit{S}}_{\mathrm{f}})$ is achieved by executing an inner product with a surface vector (S_f) on the left and right parts of Eq. (20), giving

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\mathit{\varphi }}_f^r& ={\left({\displaystyle \frac{H\left(U^m\right)}{A_P}}\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}\\ & -\left({\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right)}_f{\left(\mathit{g}\cdot \mathit{X}\right)}_f^n{\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}\mathit{\rho }\right)}^n\left|{\mathit{S}}_{\mathrm{f}}\right|\right)\\ & -{\left({\displaystyle \frac{{\left(\mathrm{\Omega }{\mathit{e}}_{\mathrm{3}}\times \mathit{U}\right)}^n}{A_P}}\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}.\end{array}\end{array}$$

Equation (21) completed the flux calculation without considering the influence of the pressure term. The pressure contribution in terms of a flux can be expressed as

$$\begin{array}{}\text{(24)}& {\left({\displaystyle \frac{-\mathrm{\nabla }{p}_{\mathit{\_}\mathrm{rgh}}}{{\mathit{\rho }}_{\mathrm{0}}{A}_P}}\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}={\left({\displaystyle \frac{-\mathrm{1}}{A_P}}\right)}_f\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}{p}_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}\right)\left|{\mathit{S}}_{\mathrm{f}}\right|.\end{array}$$

Then, Eq. (22) is added to Eq. (21) to yield

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\mathit{\varphi }}_f^{m+\mathrm{1}}& ={\left({\displaystyle \frac{H\left(U^m\right)}{A_P}}\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}\\ & -\left({\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right)}_f{\left(\mathit{g}\cdot \mathit{X}\right)}_f^n{\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}\mathit{\rho }\right)}^n\left|{\mathit{S}}_{\mathrm{f}}\right|\right)\\ & -{\left({\displaystyle \frac{{\left(\mathrm{\Omega }{\mathit{e}}_{\mathrm{3}}\times \mathit{U}\right)}^n}{A_P}}\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}\\ & -{\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right)}_f\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}{p}_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}\right)\left|{\mathit{S}}_{\mathrm{f}}\right|.\end{array}\end{array}$$

Combined with Eq. (21), Eq. (23) is simplified and rewritten as

$$\begin{array}{}\text{(26)}& {\mathit{\varphi }}_f^{m+\mathrm{1}}={\mathit{\varphi }}_f^r-{\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right)}_f\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}{p}_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}\right)\left|{\mathit{S}}_{\mathrm{f}}\right|.\end{array}$$

Using conservation of mass, we solve the pressure field $p_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}$, which results in

$$\begin{array}{}\text{(27)}& \sum _{f\in \partial V_P}{\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right)}_f\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}\mathrm{\nabla }{\displaystyle \frac{\mathrm{1}}{f}}{p}_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}\right)\left|{\mathit{S}}_{\mathrm{f}}\right|=\sum _{f\in \partial V_P}{\mathit{\varphi }}_f^r.\end{array}$$

The preconditioned conjugate gradient method is used to solve the linear system constructed by Eq. (25) (OpenFOAM, 2019). After $p_{\mathit{\_}\mathrm{rgh}}^{m+\mathrm{1}}$ is obtained using Eq. (25), we calculate the volume flux using Eq. (24) for each face. The cell-centred velocity fields ${\mathit{U}}_P^{m+\mathrm{1}}$ are calculated by reconstructing the face velocity flux using the following expression (Deshpande, 2012):

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\mathit{U}}_P^{m+\mathrm{1}}& ={\mathit{U}}_P^r+\left({\displaystyle \frac{\mathrm{1}}{A_P}}\right){\left(\sum _{f\in \partial V_P}{\displaystyle \frac{{\mathit{S}}_{\mathrm{f}}{\mathit{S}}_{\mathrm{f}}}{\left|{\mathit{S}}_{\mathrm{f}}\right|}}\right)}^{-\mathrm{1}}\\ & \cdot \left(\sum _{f\in \partial V_P}\left({\displaystyle \frac{{\mathit{\varphi }}_f^{m+\mathrm{1}}-{\left({\mathit{U}}_P^r\right)}_f\cdot {\mathit{S}}_{\mathrm{f}}}{{\left(\mathrm{1}/A_P\right)}_f}}\right){\displaystyle \frac{{\mathit{S}}_{\mathrm{f}}}{\left|{\mathit{S}}_{\mathrm{f}}\right|}}\right).\end{array}\end{array}$$

Equation (26) completes the velocity field calculation of the first iteration step in the PISO algorithm. By converting the identifier m to m+1, the next PISO iteration is completed and updates the velocity in Eq. (17) with the velocity ${\mathit{U}}_P^{m+\mathrm{1}}$ calculated from Eq. (26), thereby updating p_{_rgh}, ϕ_f and U. This procedure is performed M times to guarantee that the results of the velocity and pressure together conform to the continuity and momentum equations. Considering that PISO iteration levels require more than 1, but typically not more than 4 (OpenFOAM, 2019), we specify that the number of PISO iteration levels is 3 in the computations presented in this paper. After completing the three iterations, the converged values are considered the result of the next time step (n+1), namely,

$$\begin{array}{}\text{(29)}& {\mathit{\varphi }}_f^{n+\mathrm{1}}={\mathit{\varphi }}_f^M,{\mathit{U}}_P^{n+\mathrm{1}}={\mathit{U}}_P^M,p_{\mathit{\_}\mathrm{rgh}}^{n+\mathrm{1}}=p_{\mathit{\_}\mathrm{rgh}}^M.\end{array}$$

We discretise the convection–diffusion equation of density (Eq. 12) to obtain

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\displaystyle \frac{V_P}{\mathrm{\Delta }t}}({\mathit{\rho }}_P^{n+\mathrm{1}}-{\mathit{\rho }}_P^n)& +\sum \left({\mathit{\varphi }}_f^{n+\mathrm{1}}{\mathit{\rho }}_f^{n+\mathrm{1}}\right)\\ & =\sum k\left[\left|{\mathit{S}}_{\mathrm{f}}\right|{\displaystyle \frac{{\mathit{\rho }}_N^{n+\mathrm{1}}-{\mathit{\rho }}_P^{n+\mathrm{1}}}{\left|d\right|}}\right].\end{array}\end{array}$$

At the end of the iteration procedure, we bring the results of the volume flux into Eq. (28) to calculate the density field at the next time (${\mathit{\rho }}_P^{n+\mathrm{1}})$, thereby updating the density field for the next step calculation ($\mathrm{\Delta }t=t^{n+\mathrm{2}}-t^{n+\mathrm{1}}$).

2.3 Initialised field of ISW generation

ISW generation methods mainly include the gravity collapse mechanism, double push-pedal method (Fu et al., 2008), velocity-inlet method (Gao et al., 2012), mass source method (Wang et al., 2018), initialisation method, and methods addressing the interaction between tidal current and topography. For example, Hsieh et al. (2014) investigated the flow evolution of a depression ISW generated by the gravity collapse mechanism. Cheng et al. (2020) studied the interaction between ISWs and a cylinder using the gravity collapse mechanism. The initialisation method involves solving the internal solitary wave theory at the initial moment, such as the Korteweg–de Vries (KdV) equation (Grimshaw et al., 2010), the modified KdV (mKdV) equation, the extended KdV (eKdV) equation, the forced KdV equation, the Ostrovsky equation (Li and Farmer, 2011), the Miyata–Choi–Camassa (MCC) model (Miyata, 1985, 1988; Choi and Camassa, 1999) and the Dubreil–Jacotin–Long (DJL) equation (Long, 1953; Turkington, 1991; Brown and Christie, 1998; Dunphy et al., 2011), to obtain the wave surface and velocity field. The method of an interacting between tidal current and terrain that stimulates ISWs is adopted by many scholars, such as Farmer and Smith (1980), Lamb et al. (1994), and Shaw et al. (2009).

In this paper, the method of initialising the field is selected to generate the ISWs. To increase the application range of the ISWFoam model, two initialisation methods are provided, including solving the weakly nonlinear models of the eKdV equation (Helfrich and Melville, 2006) and the fully nonlinear models of the DJL equation for continuously stratified fluids (Turkington, 1991; Dunphy et al., 2011). The Dubreil–Jacotin–Long (DJL) equation is expressed as

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\eta }+{\displaystyle \frac{N^{\mathrm{2}}\left(z-\mathit{\eta }\right)}{c^{\mathrm{2}}}}\mathit{\eta }=\mathrm{0},& \mathit{\eta }=\mathrm{0}\mathrm{at}z=\mathrm{0},-H\\ & \mathit{\eta }=\mathrm{0}\mathrm{at}\left|x\right|\to \mathrm{\infty },\end{array}\end{array}$$

where η is the isopycnal displacement, H is the water depth, c is the propagation speed, N is the definition of the buoyancy frequency, and z is the vertical position.

$$\begin{array}{}\text{(32)}& N^{\mathrm{2}}\left(z\right)=-g{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}\left(z\right)}{\mathrm{d}z}},\end{array}$$

where ρ₀(z) is the reference density, and g is the gravitational acceleration.

By solving the above DJL equation we can obtain η and c, and then through the relationship ψ=ηc, where ψ is the stream function, we can obtain the wave-induced velocity field. We use the DJLES open-source package provided by Dunphy et al. (2011) to solve the DJL equations. Then we input the initial field calculated by DJLES into OpenFOAM to obtain the initial field required for OpenFOAM numerical simulation.

Another theory of ISWFoam model wave generation involves the weakly nonlinear models of the eKdV equation. Using the first-order stream function for the DJL equation, we can obtain the well-known KdV equation and further obtain the eKdV equation. For the specific derivation, please refer to the paper by Lamb and Yan (1996). The eKdV equation (Helfrich and Melville, 2006) is

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}+(c_{\mathrm{0}}+c_{\mathrm{1}}\mathit{\zeta }+c_{\mathrm{3}}{\mathit{\zeta }}^{\mathrm{2}}){\displaystyle \frac{\partial \mathit{\zeta }}{\partial x}}+c_{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\zeta }}{\partial x^{\mathrm{3}}}}=\mathrm{0},\text{(34)}& {\displaystyle }{c}_{\mathrm{0}}^{\mathrm{2}}={\displaystyle \frac{g{h}_{\mathrm{1}}{h}_{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{1}}\right)}{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}}},\text{(35)}& {\displaystyle }{c}_{\mathrm{1}}=-{\displaystyle \frac{\mathrm{3}{c}_{\mathrm{0}}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}^{\mathrm{2}}{h}_{\mathrm{2}}}},\text{(36)}& {\displaystyle }{c}_{\mathrm{2}}={\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{6}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{1}}^{\mathrm{2}}{h}_{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}}},\text{(37)}& {\displaystyle }{c}_{\mathrm{3}}={\displaystyle \frac{\mathrm{3}{c}_{\mathrm{0}}}{h_{\mathrm{1}}^{\mathrm{2}}{h}_{\mathrm{2}}^{\mathrm{2}}}}\left[{\displaystyle \frac{\mathrm{7}}{\mathrm{8}}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}}}\right)}^{\mathrm{2}}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}^{\mathrm{3}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}^{\mathrm{3}}}{{\mathit{\rho }}_{\mathrm{1}}{h}_{\mathrm{2}}+{\mathit{\rho }}_{\mathrm{2}}{h}_{\mathrm{1}}}}\right],\end{array}$$

where ζ is the isopycnal vertical displacement; c₀ is the linear phase speed; the coefficients c₁, c₂ and c₃ are functions of the steady background stratification and shear through the linear eigenmode (vertical structure function) of interest (Helfrich and Melville, 2006); h₁ and h₂ are the mean upper and lower layer depths, respectively; and ρ₁ and ρ₁ are the fluid densities of the upper and lower layers, respectively. The theoretical solution of Eq. (31) above is

$$\begin{array}{}\text{(38)}& {\displaystyle }\mathit{\zeta }={\displaystyle \frac{a}{B+\left(\mathrm{1}-B\right){\cosh }^{\mathrm{2}}\left[{\mathit{\lambda }}_{\mathrm{eKdV}}\left(x-c_{\mathrm{eKdV}}t\right)\right]}},\text{(39)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{eKdV}}^{\mathrm{2}}={\displaystyle \frac{a}{\mathrm{12}{c}_{\mathrm{2}}}}\left(c_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{c}_{\mathrm{3}}a\right),\text{(40)}& {\displaystyle }{c}_{\mathrm{eKdV}}=c_{\mathrm{0}}+{\displaystyle \frac{a}{\mathrm{3}}}\left(c_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{c}_{\mathrm{3}}a\right),\text{(41)}& {\displaystyle }B={\displaystyle \frac{-a{c}_{\mathrm{3}}}{\mathrm{2}{c}_{\mathrm{1}}+a{c}_{\mathrm{3}}}},\text{(42)}& {\displaystyle }{u}_{\mathrm{1}}=-c_{\mathrm{eKdV}}{\displaystyle \frac{\mathit{\zeta }}{h_{\mathrm{1}}-\mathit{\zeta }}},u_{\mathrm{2}}=c_{\mathrm{eKdV}}{\displaystyle \frac{\mathit{\zeta }}{h_{\mathrm{2}}+\mathit{\zeta }}},\end{array}$$

where a is the ISW amplitude, λ_eKdV is the wavelength, c_eKdV is the wave speed, B is an auxiliary parameter, and u₁ and u₂ are the speeds of the upper and lower layers of the fluid, respectively. The waveform and velocity field of the ISWs are solved at the initial moment by the developed function and then assigned to the calculation domain.

The vertical profile of the initial density is given by a hyperbolic tangent function profile (Aghsaee et al., 2010):

$$\begin{array}{}\text{(43)}& \stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{1}}+{\mathit{\rho }}_{\mathrm{2}}}{\mathrm{2}}}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{1}}}{\mathrm{2}}}\tanh \left({\displaystyle \frac{z-z_{\mathrm{pyc}}}{d_{\mathrm{pyc}}}}\right),\end{array}$$

where z is the vertical position; ρ₁ and ρ₂ are the fluid densities of the upper and lower layers, respectively; z_pyc is the location of the centre of the pycnocline; and d_pyc is the thickness of the pycnocline. In this paper, unless otherwise specified, the form of the density profile adopts Eq. (41). The internal solitary wave surface is obtained by calculating the gradient of the density field, and the absolute value of the maximum value of the gradient represents the vertical position of the wave surface. Notably, the density profile of the actual ocean is not always hyperbolic, so our model provides a function for users to modify the density profile according to the actual situation.

2.3.1 Comparison between the DJL equation and the eKdV equation

To compare the DJL equation and the eKdV equation, we set up a numerical simulation, which includes a tank that is 15 m long, 1 m wide and has a water depth of 0.5 m. The depths of the upper (h₁) and lower (h₂) layers are 0.1 and 0.4 m, respectively; the densities of the upper and lower layers are 1022 and 1028 kg m⁻³, respectively; the location of the centre of the pycnocline (z_pyc) is 0.4 m; the pycnocline thickness (d_pyc) is 0.04 m vertically; the initial ISW amplitude (a) is 0.065 m; and the location of the centre of ISW is 12.5 m. The ISWs propagate from right to left. The measuring point P is set at a position 10 m away from the initial ISW. The grid is uniform in the x direction, y direction and z direction, and the sizes are $\mathrm{\Delta }x=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ m, $\mathrm{\Delta }y=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ m and $\mathrm{\Delta }z=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ m, respectively. Slip boundary conditions are applied to the bottom and both sides, while cyclic boundary conditions are assigned to the inlet and outlet boundaries. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.

Figure 1 Comparison chart of the horizontal velocity component field: DJL equation (left) and eKdV equation (right).

Figure 1 shows the comparison of the horizontal velocity component field when the DJL equation and the eKdV equation are used to generate ISWs. At the initial moment, the ISW generated by the eKdV equation is not as smooth as the ISW generated by the DJL equation, and the horizontal velocity at the interface area is discontinuous as shown in Fig. 1a and b. With the propagation of ISWs, the ISWs generated by the DJL equation are always smooth at the interface area, and the velocity field is always continuous as shown in Fig. 1a, c, e and g. Correspondingly, the ISW generated by the eKdV equation gradually produces a gradient in the vertical direction of the horizontal velocity in the interface area; thus, the interface area becomes smooth, and the velocity becomes continuous. Fig. 1d shows this evolution process, which is basically completed in 5 s as shown in Fig. 1f. At 50 s, the difference between the horizontal velocity fields of the two equations is very small as shown in Fig. 1g and h.

Figure 2 shows the comparison of the vertical velocity component field when the DJL equation and the eKdV equation are used to generate ISWs. Since the theoretical solution of the eKdV equation only obtain the average horizontal velocity of the upper and lower layers of the fluid, there is no vertical velocity at the initial moment, as shown in Fig. 2b. With the propagation of ISWs, the vertical velocity field will gradually be generated and finally stabilised, and the stable time occurs at 5 s as shown in Fig. 2b, d, f and h. At 50 s, the difference between the vertical velocity fields of the two equations is very small as shown in Fig. 2g and h.

Figure 2 Comparison chart of the vertical velocity component field: DJL equation (left) and eKdV equation (right).

The ISW propagates for 10 m, and the amplitudes of the ISWs generated by the DJL equation and the eKdV equation are reduced by 9.88 % and 17.96 %, respectively, as shown in Fig. 3. Overall, the reduction in energy leads to the attenuation of the amplitude of the ISW, which in turn reduces the wave speed. Except for the difference in initial fields, the grid sizes, time step, turbulence model and other features are the same. Therefore, the initial stage of ISWs generated by the eKdV equation leads to excessive energy loss compared with those generated by the DJL equation. From the above analysis of the velocity field, we know that the method of initialising the field with the eKdV equation requires a period of movement before the jump of the velocity field develops into a field with continuous changes in velocity. In addition, the DJL equation, as a fully nonlinear model, can better reflect its superiority for internal waves with strong nonlinearity. Therefore, the wave generation of the subsequent numerical cases in this paper adopts the method of initialising the field with the DJL equation.

Figure 3 Time series of the interface displacement. The probe was 10 m away from the initial ISW.

3 Grid independence of the ISW simulation

These grid independence tests were performed in the horizontal and vertical directions by applying meshes of different sizes. The sizes of the mesh determined in this paper are calculated based on the amplitude of the ISW and a characteristic length determined through the integration of the wave profile (Michallet and Ivey, 1999):

$$\begin{array}{}\text{(44)}& L={\displaystyle \frac{\mathrm{1}}{a}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\mathit{\zeta }\left(x\right)\mathrm{d}x,\end{array}$$

where ζ is the isopycnal vertical displacement and a is the ISW amplitude.

To determine the appropriate mesh size, the propagation of ISWs on flat bottoms is calculated, and the numerical results are compared with the DJL theoretical solution. We set up a numerical simulation, which includes a tank that is 50 m long, 0.5 m wide and has a water depth of 0.5 m. The depths of the upper (h₁) and lower (h₂) layers are 0.1 m and 0.4 m, respectively; the densities of the upper and lower layers are 1000 kg m⁻³ and 1030 kg m⁻³, respectively; the location of the centre of the pycnocline (z_pyc) is 0.4 m; the pycnocline thickness (d_pyc) is 0.05 m vertically; and the ISW amplitude a is 0.065 m. The measuring point P is set at a position 10L away from the initial ISW. The sponge layer on both sides, whose length is the double wave characteristic length, has been checked to properly dissipate the reflected wave. Slip boundary conditions are applied to the bottom and both sides, while cyclic boundary conditions are assigned to the inlet and outlet boundaries. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.

3.1 Grid independence in the horizontal direction

First, we analyse the grid independence in the horizontal direction, with a constant cell height of $\mathrm{\Delta }z=a/\mathrm{20}$ m . Figure 4 shows the results of the comparison of the waveform and decay rate in the horizontal direction at probe P1 with the ISWFoam using a wide range of grid configurations. The results show a negligible difference in the waveform when the mesh size is less than $L/$40, so it is difficult to accurately analyse the grid independence just by the waveform. A traditional decay rate parameter is adopted, namely δ=(a_probe – $a_{\mathrm{initial}})/a_{\mathrm{initial}}$, where a_initial is the ISW amplitude value at the initial moment, and a_probe is the ISW amplitude value of the probe 10L away from the initial ISW. Figure 4b shows the relationship between the decay rate of the ISW amplitude and the grid quantity per unit length for different mesh sizes. As shown in Fig. 4b, the decay rate of the ISW amplitude tends to be smooth as the grid number per unit length increases to 160 ($\mathrm{\Delta }x=L/$150), and then the increase in the grid quantity has a relatively small effect on the decay rate. Therefore, for ISWFoam developed in this paper, we suggest that the dimensions of the horizontal grid are $L/$150.

Figure 4 Grid independence in the horizontal direction at probe P1: (a) comparison of waveform and (b) decay rate.

3.2 Grid independence in the vertical direction

Second, we analyse the grid independence in the vertical direction, with a constant cell width of $\mathrm{\Delta }x=L/$150 m. Figure 5 shows the results of a comparison of the waveform and decay rate of the ISW amplitude in the vertical direction at probe P1 with the ISWFoam using a wide range of grid configurations. The results also show a negligible difference in the waveform when the mesh size is less than $a/$10, so it is difficult to accurately analyse the grid independence just by the waveform. As shown in Fig. 5b, the decay rate of the ISW amplitude decreases as the grid quantity increases in a wave height range before the numerical oscillation occurs. Here, we assume that the grid size with the decay rate of the ISW amplitude less than 1 % is the appropriate vertical grid size; namely, the vertical grid size is $a/$25 m. Therefore, for ISWFoam developed in this paper, we suggest that the dimensions of the vertical grid be $a/$25.

Figure 5 Grid independence analysis in the vertical direction at probe P1: (a) comparison of waveform and (b) decay rate.

4 Model verification and results

To verify the numerical model of the ISWs, the propagation of ISWs on a flat bottom section, submerged triangular ridge and slopes is calculated, and the numerical results are compared with the corresponding experimental results. To verify the correctness of Coriolis code implantation and reflect the role of local mesh refinement, the propagation of ISWs on a flat bottom section of actual ocean scale and a Gaussian ridge is calculated.

Figure 6 Schematic diagram of probe position (P1–P5) (Hsieh et al., 2014).

Figure 7 Density contours at different moments.

4.1 ISW propagating on a flat bottom section

4.1.1 Experimental data used

In this section, ISWFoam is verified by employing ISWs propagating on a flat bottom section with Case Flat_4 in the continuously stratified laboratory experiment described in Hsieh et al. (2014). The physical dimensions and ultrasonic probe locations in the experiments of Hsieh et al. (2014), as shown in Fig. 6, are adopted to establish the numerical computation domain. We set up a numerical tank of the experiment of Hsieh and co-authors, which includes a tank that is 15 m long, 0.5 m wide and has a stable water depth of 0.5 m. The fluid densities of the upper (ρ₁) and lower (ρ₂) layers are 996 and 1030 kg m⁻³, respectively; the ISW amplitude a is 0.068 m; the location of the centre of the pycnocline (z_pyc) is 0.4 m; the pycnocline thickness (d_pyc) is 0.04 m vertically; and the depths of the upper (h₁) and lower (h₂) layers are 0.1 m and 0.4 m, respectively. The grid is uniform in the x direction, y direction and z direction, and the sizes are $\mathrm{\Delta }x=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{2}}$ m, $\mathrm{\Delta }y=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{2}}$ m and $\mathrm{\Delta }z=\mathrm{2.72}\times {\mathrm{10}}^{-\mathrm{3}}$ m, respectively. The sponge layer on both sides, whose length is double the wave characteristic length, has been checked to properly dissipate the reflected wave. Slip boundary conditions are applied to the bottom and both sides, while cyclic boundary conditions are assigned to the inlet and outlet boundaries. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.

4.1.2 Comparisons between the numerical and experimental results

Figure 7 shows the density contours at three different times from Case Flat_4 in the laboratory experiment of Hsieh and coworkers, showing the stable evolution of an ISW. The results also show the realistic evolution of the thickening of the pycnocline after ISW propagation because of convection and diffusion. At the same time, the propagation of the ISW is stable and unbroken.

To further verify the model, the waveform is compared between the numerical simulations and the experimental measurements, and the measurement point selection is the same as the experimental setting, as shown in Fig. 6. Figure 8 shows the comparison results between the waveform simulated by ISWFoam and the experimental results at probes P1–P5. Figure 8 shows that the results of the numerical simulations agree with the experimental results (red circle). Notably, the laboratory wave height at the probe P1 measurement point is greater than the numerical simulation results, and the wave surface of the laboratory wave is not smooth, which is caused by the wave generation method using the gravity collapse mechanism in the laboratory. In general, the model developed in this paper can simulate the generation and evolution of ISWs in continuously stratified fluids.

Figure 8 Comparison of the waveform between the experimental results and numerical simulation results at probes P1–P5.

4.2 ISW propagating over a submerged triangular ridge

4.2.1 Experimental data used

In this section, the validation of the numerical model is conducted through an ISW propagating over a submerged triangular ridge with the continuously stratified experiments described in Hsieh et al. (2015). The laboratory tank is 12 m long and has a stable water depth of 0.5 m, with which the fluid system has a finite thickness of the pycnocline. The specific experimental parameters used for validation of ISWFoam include the various depths of the upper (h₁) and lower (h₂) layers; the fluid density of the upper (ρ₁) and lower (ρ₂) layers of 996 kg m⁻³ and 1030 kg m⁻³, respectively; the ISW amplitude (α=0.056 m); the location of the centre of the pycnocline (z_pyc=0.4 m); the thickness of the pycnocline (d_pyc=0.04 m vertically); the height of the isosceles triangular ridge (h_s=0.30 m vertically); and the slope angle of the ridge for α=45^∘. The physical dimensions, and ultrasonic probe locations in the experiments of Hsieh et al. (2015), as shown in Fig. 9, are adopted to establish the numerical computation domain.

Figure 9 Schematic illustration of the laboratory set-up and the locations of the probes (Hsieh et al., 2015).

4.2.2 Numerical implementation

The numerical tank is designed to reproduce the experiment described in Fig. 9. The unstructured grid and local mesh refinement based on the finite-volume method are used to construct the computational domain and discretise the governing equations. The grid is uniform in the x direction, y direction and z direction, and the sizes are $\mathrm{\Delta }x=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m, $\mathrm{\Delta }y=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m and $\mathrm{\Delta }z=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ m, respectively. The precise grid described triangular ridge section is $\mathrm{\Delta }x=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ m, $\mathrm{\Delta }y=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ m and $\mathrm{\Delta }z=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ m at the slope, as shown in Fig. 10. The sponge layer on both sides, whose length is double the wave characteristic length defined through integration of the wave profile in Sect. 3 for this case, has been checked to absorb the reflected wave well. Rigid wall conditions are applied to both sides, while the slip and slip conditions are assigned to the bottom and the surface of the submerged ridge boundaries, respectively. The top boundary is a rigid lid. The inlet and outlet boundaries adopt cyclic boundary conditions. The boundary conditions related to the density field are no-flux boundary conditions.

Figure 10 Schematic of the mesh.

Figure 11 Comparison of the waveform between the experimental results in Hsieh et al. (2015) and numerical simulation results at probes P1–P5.

4.2.3 Comparisons between the numerical and experimental results

Figure 11 shows the comparison results between the waveform calculated by ISWFoam and the experimental results at probes P1–P5. In each subplot, the results of the numerical simulations (blue line) are found to be in good agreement with the experimental results (red circle). From Fig. 11a, the numerical simulation result of the probe P1 measurement after 20 s is different from the experimental results, which is caused by the different ISW generation methods. For the experimental results, the first large leading ISW is formed via the gravity collapse mechanism, which is trailed by a train of small-amplitude mode-one waves that is generated due to shear instabilities. However, the initialisation method used to generate an initial ISW for the numerical simulation in this paper is more stable than the gravity collapse mechanism, so the rear part of the ISW is relatively flat compared to the experimental results for probe P1. In Fig. 11, the waveform of the ISW gradually evolves towards a flat waveform due to the interaction between the ISW and the ridge. In general, the model developed in this paper can simulate the interaction between ISWs and structures.

4.3 ISW propagating on a slope

To verify the ability and accuracy of simulating the ISW breaking of the numerical model, two continuously stratified laboratory experiments (12 and 15) of Michallet and Ivey (1999) are chosen for the simulation in this section. The experimental set-up is represented schematically in Fig. 12. We set up a numerical tank of the experiment of Michallet and Ivey (1999), which includes a tank that is L=4.2 m long, 0.25 m wide and has a water depth of 0.15 m. The layer thickness ratio $(h/H)$ varies from 0.60 to 0.91. A linear slope s=0.214 starts at 0.7 m from the right end of the tank for experiments 12 and 15. The grid is uniform in the x direction, y direction and z direction, and the sizes are $\mathrm{\Delta }x=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ m, $\mathrm{\Delta }y=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$ m and $\mathrm{\Delta }z=\mathrm{1.25}\times {\mathrm{10}}^{-\mathrm{3}}$ m, respectively. The precise grid describing the slope section is $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{6.25}\times {\mathrm{10}}^{-\mathrm{4}}$ m and $\mathrm{\Delta }z=\mathrm{3.125}\times {\mathrm{10}}^{-\mathrm{4}}$ m at the slope.

Figure 12 Schematic diagram of the laboratory set-up.“C” and “US” represent the experimental devices at the probes.

Figure 13 Time series of the interface displacement. The probe was 99.8 cm away from the start of the slope.

The sponge layer on the left side, whose length is double the wave characteristic length, is checked to properly dissipate the reflected wave. Slip boundary conditions are applied to the bottom and both sides, while slip boundary conditions are assigned to the top boundaries. The boundary conditions related to the density field are no-flux boundary conditions.

The vertical profile of the initial density is given by a hyperbolic tangent function profile

$$\begin{array}{}\text{(45)}& \stackrel{\mathrm{‾}}{\mathit{\rho }}\left(z\right)={\mathit{\rho }}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{\mathrm{2}}}\left\{\mathrm{1}+\tanh \left[{\displaystyle \frac{-\left(z-z_{\mathrm{pyc}}\right)}{d_{\mathrm{pyc}}}}\right]\right\},\end{array}$$

where z is the vertical position, ${\mathit{\rho }}_{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{3}}$ kg m⁻³ is the base density field, Δρ is the change in the density, z_pyc is the location of the centre of the pycnocline, and d_pyc is the thickness of the pycnocline.

4.3.1 Case one and results

The first case of model verification is experiment 12 of Michallet and Ivey (1999) in this section. The layer thickness ratio $(h/H)$ is 0.84, and the density change (Δρ) value is 14 kg m⁻³. Figure 13 presents the time series for the interface displacement (ζ) for experiment 12 of Michallet and Ivey (1999). The results indicate reasonably good agreement between the time series of the simulated interface displacement and that of the laboratory results. The first trough centred around t=25 s represents the incident ISW propagating the probe 99.8 cm away from the start of the slope. The second trough centred at approximately t=87 s represents the reflected ISW at the generation side, which has a smaller amplitude and a longer wavelength than the incident ISW as the energy in the wave decreases. As shown in Fig. 13, the smooth waveform of the incident ISW of the numerical simulation indicates that the initialisation method of wave generation in this paper is more stable than the experiment.

Figure 14 Comparison of the velocity field between the experimental observation results in Michallet and Ivey (1999) (left) and numerical simulation results (right).

Figure 14 shows a comparison of ISWFoam results and the experimental observations of the velocity field associated with the ISW run-up process along the slope. The model effectively reproduces laboratory tests, such as the intensity and direction of the velocity field, the location of the vortices and the occurrence of boundary-layer separation beneath the ISW. Therefore, the model developed in this paper can reflect the ISW breaking phenomenon during the propagation of ISWs along the slope.

4.3.2 Case two and results

Another laboratory experiment that more clearly shows the ISW breaking phenomenon from Experiment 15 of Michallet and Ivey (1999) is used to verify the numerical model presented in this paper, and the corresponding numerical case is set corresponding to it. The layer thickness ratio $(h/H)$ is 0.77, and the density change (Δρ) value is 47 kg m⁻³. The wave amplitude and phase velocity at the slope calculated by ISWFoam are $\mathrm{2.71}\times {\mathrm{10}}^{-\mathrm{3}}$ m and $\mathrm{10.83}\times {\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹, which fit well with the experimental results of $\mathrm{2.7}\times {\mathrm{10}}^{-\mathrm{3}}$ m and $\mathrm{10.8}\times {\mathrm{10}}^{-\mathrm{1}}$ m s⁻¹.

Figure 15 shows the results of the numerical simulations of the ISWs propagating along the slope and wave breaking using ISWFoam. As the ISW propagates to the slope, according to the conservation of mass, the upper fluid forward and the downward velocity of the lower fluid increasing along the slope results in the formation of a thin boundary layer, as shown in Fig. 15a, b and c. At the same time, the amplitude of the ISW increases, and the rear of the ISW gradually becomes very steep but does not overturn. With the development of the ISW, the rear waveform of the ISW cannot maintain its stability and overturns forward, resulting in wave breaking, as shown in Fig. 15d. After wave breaking occurs, the denser lower layer flow accelerates into the less dense upper layer flow, forming a mixture region, as shown in Fig. 15e. After the lower layer flow is drawn downward from beneath the ISW, a mixing region comprised of vortices is pushed upwards along the slope while the leading waveform is reflected, as shown in Fig. 15f, g and h. Figure 15i and j show the falling process of ISWs. From the perspective of the entire process of wave breaking, the steepening of the rear waveform in this case is the main reason for wave breaking.

Figure 15 Temporal and spatial variations in the ISWs breaking calculated using ISWFoam (the black line represents the waveform).

For comparison with the flow visualisation image of the experiment, a specified thickness of the pycnocline is presented, and the pycnocline ranges from 1003 to 1045 kg m⁻³ with dark colours as shown in Fig. 16.

Figure 16 Comparison of the density fields between the experimental observation results in Michallet and Ivey (1999) (left) and the numerical simulation results (right).

Figure 16 compares the ISWFoam results and the experimental results of Michallet and Ivey (1999) before, during and after ISW breaking. The results indicate that some main features of the laboratory tests are reasonably well reproduced by ISWFoam, such as the profile of ISW, the location of the wave breaking point, ISW arrival time, and spatial and temporal changes in the mixture region. Therefore, the model developed in this paper can accurately simulate the ISW breaking phenomenon during the propagation of ISWs along the slope.

5 Simulation examples at the field scale

The ISWFoam model developed in the present paper can be used as a tool to investigate the interaction between ISWs and complex structures and topography. In this section, two numerical examples are presented to show the capability of ISWFoam on field-scale simulation.

5.1 ISW propagating over a 3D Gaussian ridge

We designed a case of an ISW propagating over a 3D Gaussian ridge. The 3D Gaussian ridge is obtained by rotating a 2D Gaussian ridge:

$$\begin{array}{}\text{(46)}& z=a{e}^{-{\left(x/l\right)}^{\mathrm{2}}},\end{array}$$

where a is the ridge amplitude, and l is the standard deviation. With a=100 m and l=10, we can obtain a 2D Gaussian ridge with a height of 100 m and a bottom width of 40 m. Subsequently, the 3D Gaussian ridge can be obtained after a vertical rotation of 180^∘.

Figure 17 Schematic of the local refinement of the grid.

Figure 18 Temporal and spatial variation in the ISWs propagating over a 3D Gaussian ridge.

We set up a 3D numerical tank, which includes a tank that is 3 km long, 400 m wide (y direction from −200 to 200 m) and has a water depth of 120 m. The depths of the upper (h₁) and lower (h₂) layers are 20 and 100 m, respectively; the densities of the upper and lower layers are 1000 and 1030 kg m⁻³, respectively; the location of the centre of the pycnocline (z_pyc) is 100 m; the pycnocline thickness (d_pyc) is 1.5 m vertically; and the ISW amplitude a is 20 m. The Gaussian ridge is located at 800 m horizontally. The grid is gradually changed from Δx=20 m to Δx=2.5 m in the x direction, the grids in the y direction are uniform with a constant cell width of Δy=2.5 m, and the grids in the z direction are non-uniform, with a minimum cell height of Δz=1 m near the interface of the ISW. The precise grid described the 3D Gaussian ridge section as $\mathrm{\Delta }x=\mathrm{3.9}\times {\mathrm{10}}^{-\mathrm{2}}$ m, $\mathrm{\Delta }y=\mathrm{3.9}\times {\mathrm{10}}^{-\mathrm{2}}$ m and $\mathrm{\Delta }z=\mathrm{1.56}\times {\mathrm{10}}^{-\mathrm{2}}$ m, as shown in Fig. 17. The sponge layer on both sides, whose length is the double wave characteristic length, has been checked to properly dissipate the reflected wave. Slip boundary conditions are applied to the bottom and both sides, while cyclic boundary conditions are assigned to the inlet and outlet boundaries. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.

Figure 18 shows the temporal and spatial variations in the ISWs propagating over a 3D Gaussian ridge. The ISW reaches the Gaussian ridge, causing the wave surface in front of the ridge to decrease, and the wave surface behind the ridge to climb up the ridge, as shown in Fig. 18a. Due to being obstructed by the Gaussian ridge, flow around a ridge and wave surface uplift are generated on both sides of the Gaussian ridge (perpendicular to the direction of wave propagation), as shown in Fig. 18b. As the ISW propagated over the Gaussian ridge, the wave surface climbed along the ridge, and at the same time, low velocity was generated behind the ridge, as shown in Fig. 18c. Since the top of the ridge is in the pycnocline, there will be a low-velocity area behind the ridge for a period of time after the ISW passes, as shown in Fig. 18d. In general, the ISWFoam model with unstructured grids and local mesh refinement can simulate the interaction between ISWs and complex structures and topography at the field scale.

5.2 ISW propagating over a hyperbolic tangent terrain

The propagation of ISWs to the shore is bound to be affected by the continental shelf, and shallow water evolution phenomena such as nonlinear evolution, breaking phenomena and waveform inversion occur on the undulating continental shelf. For simplicity, this section simplifies the continental shelf into a hyperbolic tangent terrain, and the terrain profile formula is as follows (Lamb, 2002):

$$\begin{array}{}\text{(47)}& z={\displaystyle \frac{a_T}{\mathrm{2}}}\left(\mathrm{1}+\tanh \left({\displaystyle \frac{x}{L}}\right),\right)\end{array}$$

where a_T is the ridge amplitude, which is 60 m; l is the width of the continental shelf change area, which is 200 m; and x and z are the horizontal and vertical coordinate positions, respectively.

Figure 19 Schematic of the local refinement of the grid.

Figure 20 Velocity field diagram of ISW propagation (the black solid line is the iso-density contours, t=800 to 1300 s).

We set up a 3D numerical tank, which includes a tank that is 7.5 km long, 200 m wide (y direction from −100 to 100 m) and has a water depth of 100 m (z direction from 0 to 100 m). The depths of the upper (h₁) and lower (h₂) layers are 20 and 80 m in the deep water, respectively, and the depths of the upper (h₁) and lower (h₂) layers are both 20 m in the shallow water. The densities of the upper and lower layers are 1000 and 1012 kg m⁻³, respectively; the pycnocline thickness (d_pyc) is 2 m; the location of the pycnocline (z_pyc) centre is at 80 m; and the ISW amplitude a is 20 m. The starting point of the hyperbolic tangent terrain is located at 0 m horizontally, and the top of the terrain is 40 m underwater. Therefore, ISWs will gradually propagate from a deep water area with a water depth of 100 m to a shallow water area with a water depth of 40 m. The horizontal grid is gradually changed from 40 m at the inlet to 4 m in the study area, and then from the study area to 60 m at the outlet, the grid is uniform in the y direction and z direction, which are 4 and 1 m. The hyperbolic tangent terrain is characterised by 4-level local mesh refinement, and the minimum grid size is 5 cm, as shown in Fig. 19. The sponge layer on both sides, whose length is double the characteristic wave length, has been checked to properly dissipate the reflected wave. Cyclic boundary conditions are assigned to the inlet, outlet and both sides, while the slip and non-slip conditions are assigned to the bottom and the shelf topography. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.

Figure 21 Velocity field diagram of ISW propagation (the black solid line is the iso-density contours, t=1400 to 3500 s).

Figure 22 Vorticity field (areas marked A–E represent transverse sections).

Figure 20 shows the waveform and velocity field when the ISW passes through the hyperbolic tangent terrain. From Fig. 20, it can be seen that the ISW breaks and has a significant waveform inversion when propagating from deep water of 100 m to shallow water of 40 m. As the ISWs propagate to the continental shelf, the water depth gradually becomes shallower, and the thickness of the lower fluid gradually decreases as shown in Fig. 20a. Due to the presence of the continental slope, the nonlinearity of the ISWs becomes stronger, and the trough velocity of the ISWs is significantly lower, which causes the waveform at the rear of the ISW to become steep, as shown in Fig. 20b. At the same time, the front waveform of the ISW gradually becomes flat and parallel to the shelf topography. As the waveform at the rear of the ISW becomes steeper and loses balance, the waveform at the rear of the ISW rolls forward, leading to the occurrence of ISW breaking phenomena, as shown in Fig. 20c. It is worth noting that the ISW breaking occurs at the rear of the ISW, while the front waveform does not break but transforms into another form of wave (referred to as the head wave) and continues to propagate steadily along the continental shelf. The breaking of ISWs causes severe disturbances in the water and excites a series of secondary waves at the tail of the head wave, represented by elevation ISWs (as shown in Fig. 20d, e), and then elevation ISWs propagate forward steadily in shallow water (as shown in Fig. 21g, h, i). The shelf slope of the case in this section is the same as the shelf slope of s8_c1c case studied by Lamb and Xiao (2014), both of which are 0.1. The research results of Lamb and Xiao (2014) show that waveform inversion of a depression ISW will occur at this shelf slope, and a series of elevation ISWs will be generated and propagate stably in shallow water. The simulated results have good agreement with that of Lamb and Xiao (2014).

Figure 23 Vorticity field of transverse sections.

The vortex structure has an important influence on the material transport at the bottom of the shelf, so it is very necessary to study the vortex structure when the ISW breaks. Figure 22 show the vorticity field of the ISW at the breaking stage. With the occurrence of ISW breaking, a significant anticlockwise vortex structure is generated below the waveform at the rear of the ISW, as shown in Fig. 22a. With the propagation of the head wave, the vortex climbs along the shelf, the vortex continues to develop horizontally and vertically during the upward climb, and the vertical scale is about one-third of the local water depth (as shown in Fig. 22b, c). As the vortex structure continues to climb, the vorticity decays, and the vortex structure gradually disappears, as shown in Fig. 22d and e. Combined with the velocity field in Fig. 20, it can be seen that the vorticity before and after the ISW breaks is the largest, and the vortex structure is the most obvious. As the wave train of elevation ISWs propagates steadily, the vortex structure climbs up the shelf and gradually disappears.

Figure 24 Velocity vector field of wave train in shallow water areas.

Figure 25 Vorticity field of wave train in shallow water (the black solid line is the waveform).

In order to analyse the vortex of the three-dimensional structure, Fig. 23 show the vorticity field diagram on the transverse section, and the position of the transverse section corresponds to the marked section in Fig. 22. It can be seen from Fig. 23 that the vortex also evolves in the transverse section. Obviously, the bottom vortex structure generated by ISW breaking shows three-dimensional non-uniform features.

The velocity vector field of the head wave and the wave train of elevation ISWs in the shallow water area are shown in Fig. 24. The head wave generated by the breaking of the ISW loses the original wave shape of the ISW, the wave height becomes smaller, the wavelength becomes longer, and the velocity field is still in the form of upper layer forward and lower layer backward (as shown in Fig. 24a, b). In Fig. 24c and d, the velocity field and waveform of the entire wave train following the head wave are stable, the velocity field of each wave is backwards in the upper layer and forwards in the lower layer, and the wavelength gradually becomes longer as the wave train propagates. As the wave train propagates in shallow water, there is a large vorticity in the crest and trough areas of each wave, and it propagates forward steadily as the wave train propagates, as shown in Fig. 25. Generally, the waveform inversion and breaking phenomenon of ISWs is well indicated, and the propagation and evolution of the wave train generated by waveform inversion is also accurately described through ISWFoam simulation.

6 Conclusions

In this paper, a numerical model referred to as ISWFoam with a modified k–ω SST model, established by combining the density transport equation with a fully three-dimensional (3D) Navier–Stokes equation, is developed to simulate ISWs in continuously stratified, incompressible, viscous fluids based on the finite-volume method with unstructured grids and local mesh refinement of OpenFOAM. ISWFoam provides two initial wave generation methods to generate an ISW in continuously stratified fluids, including solving the weakly nonlinear models of the eKdV equation and the fully nonlinear models of the DJL equation. The verification process presents several applications, such as ISWs propagating on flat bottoms including laboratory scale and actual ocean scale, and ISWs over a submerged triangular ridge, a Gaussian ridge and slopes. The following conclusions were obtained as a result of this study.

ISWFoam using the finite-volume method with unstructured grids and local mesh refinement can accurately simulate the generation and evolution of ISWs, the ISW breaking phenomenon, waveform inversion of ISWs and the interaction between ISWs and complex structures and topography. The method of initialising the ISW using weakly nonlinear eKdV equation models requires a period of movement before the jump of the velocity field develops into a field with continuous changes in velocity. The DJL equation wave generation method that considers the vertical velocity and the horizontal velocity along the vertical gradient is better than the eKdV equation wave generation method that only provides the horizontal average velocity. Using ISWFoam to simulate an ISW with infinite wave length, the metric for the appropriate mesh size is given as follows: the dimensions of the horizontal grid are 1$/$150th of the characteristic length, while the vertical grid takes 1$/$25th of the ISW amplitude.