Two-Layer Non-Hydrostatic Model for Generation and Propagation of Interfacial Waves

When pycnocline thickness of ocean density is relatively small, density stratification can be well represented as a two-layer system. In this article, a depth integrated model of the two-layer fluid with constant density is considered, and a variant of the edge-based non-hydrostatic numerical scheme is formulated. The resulting scheme is very efficient since it resolves the vertical fluid depth only in two layers. Despite using just two layers, the numerical dispersion is shown to agree with the analytical dispersion curves over a wide range of kd, where k is the wave number and d the water depth. The scheme was tested by simulating an interfacial solitary wave propagating over a flat bottom, as well as over a bottom step. On a laboratory scale, the formation of an interfacial wave is simulated, which also shows the interaction of wave with a triangular bathymetry. Then, a case study using the Lombok Strait topography is discussed, and the results show the development of an interfacial wave due to a strong current passing through a sill.


Introduction
Internal waves are phenomena commonly observed in the ocean. Even though they occur within the ocean body, their presence can be visible from satellite images as stripes of rough surfaces. As reported (Maxworthy, 1979;Osborne and Burch, 1980;Helfrich and Melville, 1986;Lamb, 1994), internal waves play an important role in the transport of momentum and energy within the ocean. Also, they have profound effects on the Earth's climate and ocean ecosystems.
Stratified ocean density supports the existence of internal waves within the fluid medium. Typical ocean density is low at the surface and increases as depth increases until it reaches a pycnocline layer. From there onwards ocean density is relatively constant. When pycnocline thickness is relatively small, internal waves propagate horizontally along the pycnocline, and they are called interfacial waves. In this case, ocean density stratification can be well represented as a two-layer fluid system.
In modelling wave transformation from deep water to intermediate water, dispersion effect plays an important role. Dispersion means that waves of different wavelengths travel at different phase speeds. Basically, there are two ap-proaches for dispersive wave models: the Boussinesq-type models and the non-hydrostatic models. Numerous Boussinesq models with different accuracies can be found (Nwogu, 1993;Liu et al., 2018;Lawrence et al., 2018) and the references therein. These dispersive models incorporate the non-hydrostatic pressure with the aim of having a wide range of applicability (in terms of kd, where k is the wave number and d the water depth), often with the expense of simplicity and efficiency. In this article, we focus on the non-hydrostatic models which decompose the pressure into the hydrostatic pressure and the non-hydrostatic pressure. This method has been successfully used in many models such as MITgcm (Marshall et al., 1997), non-hydrostatic ROMS (Kanarska et al., 2007), Delft3D (Bijvelds, 2003), and SWASH (Zijlema et al, 2011), see also (Fang et al, 2015;Cheung, 2012, 2018;Yamazaki et al., 2008). Among these non-hydrostatic models, the most efficient and accurate numerical method is the one proposed by Stelling and Zijlema (2003). The main breakthrough was the use of the Kellerbox scheme, which is an edge-based method, so that zero pressure along the free surface can be imposed exactly and good dispersion can be achieved by applying only two or three layers on the vertical axis. Commonly used methods are cell centered that need more than 10 layers, which is expensive in computations (Casulli and Stelling, 1998;Zhou and Stansby, 1999). Further discussions on the scheme can be found (Pudjaprasetya et al., 2017;Tjandra and Pudjaprasetya, 2015;Pudjaprasetya and Tjandra, 2014;Magdalena, 2015). In this paper, the efficient and accurate two-layer non-hydrostatic scheme is implemented to study interfacial waves, focusing on their interaction with current and topography, as well as on the generation process.
As mentioned before, our approach uses a two-layer fluid system to represent ocean density stratification, and we start with the two-layer hydrostatic model which is formulated in Section 2. Next, the hydrodynamic pressure is incorporated to yield the two-layer non-hydrostatic model. This model can directly be solved using the staggered scheme. Dispersion relation of this two-layer nonhydrostatic scheme is shown to conform to the analytical dispersion. All these aspects are described in Section 3. Validity of the scheme is discussed in Section 4, via simulations of interfacial wave dynamics: a solitary wave propagates undisturbed in shape with constant velocity, and a solitary wave propagates over a bottom step. A laboratory scale simulation of an internal wave which is generated, propagates, and interacts with a triangular barrier is performed, and finally the development of interfacial wave as a tidal current interacts with Lombok Strait topography is presented in Section 5.

Two-layer hydrostatic model
When pycnocline thickness is relatively small, ocean density can be represented reasonably well as a two-layer fluid system. Here, we consider a two-layer of homogeneous fluid, bounded above with a free surface and below with an impermeable bottom, a sketch of two-layer fluid system is depicted in Fig. 1. Under the assumption of long waves in a relatively shallow region, depth integrated twolayer formulation is suitable, which is as follows: (1) where η 1 and η 2 denote the deviation of free surface and in-terface, respectively. Depth-averaged horizontal velocities for the upper and lower layer are denoted as u 1 and u 2 , respectively, whereas fluid densities are denoted as ρ 1 and ρ 2 . All dependent variables η 1 , η 2 , u 1 , and u 2 are functions of spatial variable x and time t. Eqs. (1)-(4) are the governing equations for the two-layer hydrostatic fluid system. In the following, we formulate the dispersion relation of the linearized form of Eqs. (1)-(4) as follows: where d 1 and d 2 are the undisturbed thicknesses of the upper and lower layers, respectively. Each dependent variable is assumed to have an individual wave like form with frequency ω and wave number k as follows: where A m and U m represent the amplitudes of η m and u m , respectively, for m = 1, 2. Substituting the above ansatz into Eqs. (5)-(8) will give us a system of linear equations with four unknowns A m , U m , for m = 1, 2. After some algebra we can simplify the matrix coefficient as follows: . Dispersion relation of this two-layer hydrostatic system can be obtained by equating det(M) = 0 or Two solutions for ω/k of Eq. (10) are ⩽ with d = d 1 + d 2 . This system has two intrinsic phase velocities, c 0 corresponds to the fast mode and c 1 corresponds to the slow mode. In most cases c 0 > c 1 , hence stability condition for the hydrostatic two-layer model reads as c 0 Δt/Δx 1.

Two-layer non-hydrostatic model
The non-hydrostatic approach of Stelling and Zijlema (2003) is applied to our two-layer system, and the formulation is discussed in this section. When the hydrodynamic pressure P is incorporated, the depth integrated momentum Eqs. (3) and (4) become: In the formulations above P 1 , P 2 , and P 3 denote the hydrodynamic pressures along the surface, the interface, and along the bottom, respectively. Adopting the Keller-box approximation for pressure derivative on the right hand side of Eqs. (13) and (14), is actually to assume a linear variation of P within each layer. The terms on the right hand side of Eqs. (13) and (14) are the approximation of the depth averaged pressure derivative over the upper and lower layers respectively. In the above formulation, hydrodynamic pressure in the upper layer is approximated by (P 1 + P 2 )/2, and the same also applies to the lower layer. Applying the Keller-box scheme for the momentum-z equation for each layer, and after neglecting the non-linear terms, reads Further, incompressibility for each layer yields Further simplification is obtained if we take w 3 = 0 which comes from impermeability of the bottom boundary, and P 1 = 0 which is because hydrodynamic pressure along the free surface is zero. Hence, the two-layer non-hydrostatic model is Eqs. (1)-(2) and Eqs. (13)-(18), with eight computed variables η m , u m , w m , P m+1 , with m = 1, 2. Fig. 2 describes the locations of the computed variables. In this arrangement the grid points for w m , m = 1, 2, are located along the surface and interface, whereas grid points for P k+1 , k = 1, 2, are located along the interface and bottom topography. Different from commonly non-hydrostatic models which are cell centered, in our model variables w m and P m+1 , m = 1, 2 are located on the edge of Arakawa grid cells, which allows us to systematically set the vanishing pressure along the surface P 1 = 0, as well as the vanishing vertical velocity along the bottom w 3 = 0, which makes this model very efficient. As discussed by Stelling and Zijlema (2003), this edge based scheme allows a very small number of layers (in the order of 1-3 layers) for the accurate description of relatively short waves. Whereas the cell centered methods would typically require more than 10 layers.
⩽ ⩽ For a computational domain 0 x L, application of the staggered grid partition is as follows In this staggered scheme u 1 , u 2 are calculated at half points, whereas η 1 , η 2 , (hence h 1 , h 2 ) are calculated at full points, see Fig. 2. For time integration, a stable scheme will be achieved if the mass conservation is integrated by forward time, and the momentum-x equations are backward time, or the other way around. The full nonhydrostatic discrete system that we implemented here is as follows: Fig. 2. Staggered grid and positions of the unknowns. Note that each layer thickness h m can be computed from for m=1, 2.
(26) In the above formulations, notation * h m means upwind approximation depending on u m , and that holds for m = 1, 2. Moreover, the advection terms adv m, j+1/2 , for m = 1, 2 are approximated using momentum conservative scheme as recorded in Stelling and Duijnmejer (2003 3. Corrector: To show the differences between computations using non-hydrostatic and hydrostatic model, in Section 4, we will compare the non-hydrostatic results with those of hydrostatic. And the hydrostatic numerical scheme is just Eqs. (19) and (20) for mass conservation, and Eqs. (21) and (22) with zero right hand sides, for momentum conservation.

Dispersion relation
For the rest of this paper, we will use the non-hydrostatic scheme. So here, we first discuss the dispersion relation of the linearized form of the discrete model (19) into Eqs. (19)-(26) to yield a system of linear equations with eight unknowns A m , U m , W m , P m+1 , for m = 1, 2. After some algebra we obtain the matrix coefficient below and of the relations Eq. (35), a dispersion relation of the two-layer non-hydrostatic scheme can be obtained as roots of Eq. (35) with the help of algebraic computer software. Formulae are lengthy and tedious, therefore they are not given here. For a specific case when ρ 2 = ρ 1 = ρ, and d 2 = d 1 = d/2, Eq. (35) is reduced to the following dispersion relation of the two-layer model of Stelling and Zijlema (2003) For a specific choice of ε= 0.01, and d 1 = d/8, d 2 = 7d/8, with d = 1, in Fig. 3 the numerical dispersion curves are plotted together with the analytical dispersion curves. Note that the exact dispersion is given by (36) Just like hydrostatic model, this non-hydrostatic model also has two modes: fast mode and slow mode. The fast and slow modes correspond to the formula for ω 2 /(gk) as in Eq. (37) with + and -algebraic operations, respectively. Fig. 3 shows that the dispersion relation of our discrete model nicely approximates the exact dispersion relation (37) for a wide range of kd. Moreover, it introduces an anomaly in the dispersion relation around kd = 2.5: overestimation of the phase speed in the shallow and intermediate depth range, and underestimation in deep water. This feature is similar to the dispersion relation of stratified model (Bai and Cheung, 2012). Furthermore, for parameters chosen here, the numerical dispersion relation (35) approximates (37) for 0 < kd < 10 with an error smaller than 4%.

Numerical validation
In this section, we conduct numerical validation of the two-layer non-hydrostatic scheme Eqs. (19)-(26). We start with a simulation of a baroclinic interfacial wave.

Simulation of baroclinicity
In this subsection, we conduct a baroclinic mode oscillation on a computational domain [0, L = 32 m]. We take the following two-layer configuration: the upper layer d 1 = 1 m, ρ 1 = 1000 kg/m 3 , and the lower layer d 2 = 7 m, ρ 2 = 1010 kg/m 3 . Here we take the initial interface η 2 (x, 0) = Acos(kx), with amplitude A = 0.01 m and wave number k = 4π/L and the initial surface in order to have a hot start for a baroclinic oscillation. The chosen interfacial wave is linear but dispersive, since A=d = 0.00125 and kd =π. Computations were conducted using non-hydrostatic and hydrostatic model, with Δx = 0.4 m, and Δt = 0.01 s. For both computations, wave signals at x = 0 were recorded, and the result is plotted in Fig. 4. Periods of these wave signals were measured, and the interfacial wave period we obtained is 42 s, whereas the surface wave period is 3.2 s. By using hydrostatic model, the period of wave signals is 36 s for the interfacial wave, and 1.8 s for the surface wave.
Clearly, the hydrostatic model mistakenly calculates the wave periods in this dispersive wave case. This simulation also ensures that interface oscillations are accompanied by small amplitude oscillations on the surface. This explains the rigid lid assumption that are often adopted in numerous research about internal waves (Choi and Camassa, 1999;18, 19, Grimshaw et al., 200718, 19, Grimshaw et al., , 2008. Under this rigid lid assumption, KdV-type model and their soliton solution arises (see Pudjaprasetya (2009) and the references therein). But note that we do not apply the rigid lid assumption here.

Simulation of solitary interfacial waves
Further, we conduct a simulation of solitary interfacial wave. For that purpose, the Gardner solitary wave is used as the initial hump (38) Fig. 3. Curves of numerical dispersion relation (35) in comparison with the exact dispersion relation (37). The slow mode is multiplied by ten for graphical purposes. Those curves are plotted using parameters ε= 0.01, with d 1 = d/8, d 2 = 7d/8, with d = 1.

Fig. 4.
Wave signal recorded at x = 0 in the baroclinic mode oscillation, calculated using the two-layer non-hydrostatic and hydrostatic models.
Solitary wave simulation is conducted in a two-layer system of homogeneous fluid, upper layer d 1 = 4 cm, ρ 1 = 1000 kg/m 3 , and lower layer d 2 = 28 cm, ρ 2 = 1010 kg/m 3 . Taking the parameter γ= 2.39 for the initial wave will result in a solitary wave with amplitude 1 cm. For computation we use Δx = 0.02 cm, and Δt = 0.005 s, and the results are plotted in Fig. 5. It is shown that the interfacial solitary wave propagates undisturbed in shape, with constant velocity 0.0597 m/s to be compared with c 0 = 0.0583 m/s, which is the linear phase speed of interfacial wave.
Further, we use the same solitary wave but now we let it run over a step of 10.8 cm height. Simulation results are shown in Fig. 5. Owing to a bottom step up, the solitary wave gets a ripple tail trailing backwards. In this case, no clearly reflected wave is observed, since in this case the reflection coefficient as predicted by the linear theory is , whereas the transmission coefficient is . Our finding here is in a good agreement with solitary wave simulation at a bottom step calculated using the non-hydrostatic model of the full non-linear Navier-Stokes equations in Maderich et al. (2010).

Internal wave upon a triangular barrier
In this subsection, we conduct a simulation to mimic the internal wave experiment of Sutherland et al. (2015). In a wave tank with the dimensions 197.3 cm long, 17.3 cm wide, and 40 cm high, a triangular barrier with the dimensions L b = 48 cm and the height H b = 15 cm are placed at the bottom of the tank. Salt water with the depth d 2 and density ρ 2 laid at the bottom of the tank, on the top of it is a layer of freshwater with d 1 in depth, and density ρ 1 . Initially, the left part of the tank length L l was set up so that it consists of the same two fluid layers, but the fresh fluid has depth H l , with H l > d 1 . This initial situation is depicted in Fig. 6 (top). In the experiment, this difference in height was maintained by placing a membrane.
For simulation, we use parameters d 1 = 5.5 cm, d 2 = 21.5 cm, ρ 1 = 1.000 g/cm 3 , ρ 2 = 1.050 g/cm 3 , and g = 1000 cm/s 2 . Computation is conducted using Δx = 1 cm, Δt = 0.002 s. Soon after the computation started, an internal wave in the form of a depression solitary wave is generated. It is shown that the wave interacts with the topography and breaks when it hits the submerged triangular barrier. Fig. 6 displays that our numerical simulations qualitatively resemble the experimental results of Sutherland et al. (2015). There is an obvious difference for the simulation at time t = 12 s, i.e. when the wave breaks, one reason is that our two-  layer scheme is a depth-integrated model that cannot describe breaking wave. Another reason is that, some discrepancy may be due to the inaccuracy of taking movie-snapshot for the selected time.

Onset of internal waves in Lombok Strait
As recorded in Susanto et al. (2005), internal waves were observed in Lombok Strait. Internal waves propagating northward and southward were seen from satellites images as an alternating series of bright and dark bands. The width and length of the wave-packet are approximately 85 km and 100 km, respectively. The leading wave has a wavelength of 5 km, and the rear wave has a wavelength of 2 km.
Lombok Strait is one of the outflow straits of the Indonesian Troughflow. Through Lombok Strait, warm water from the Pacific Ocean is transported to the Indian Ocean. As illustrated in Fig. 7, further south of Lombok Strait, this flux of water undergoes a narrowing passage because of Nusa Penida Island. On the right channel between Nusa Penida and Lombok, there is a complex bathymetry with a sill that further reduces the cross section of the flow (see Fig. 7). This combination of stratified water, strong tidal currents, and narrowing passage provides a favorable situation for the generation of internal waves.
Here we conduct a one dimensional study of interfacial wave generation in Lombok Strait. A cross section of Lombok Strait topography with a sill is shown in Fig. 8. It is shown that Lombok Strait has a typical depth of 800-1000 m. In the upper layer of the ocean body between the surface and pycnocline, there was a strong current flowing towards the South. Under this setting, we simulate the onset of interfacial wave in Lombok Strait using the two-layer non-hydrostatic scheme.
The two-layer setup used in the computations is as follows: the upper layer thickness is d 1 = 150 m, and fluid densities are ρ 1 = 1020 kg/m 3 , and ρ 2 = 1030 kg/m 3 for the upper and lower layer, respectively. As the generation mechanism, a horizontal current is introduced in the upper layer u 1 (x, 0) = -0.2 m/s, but no current in the lower layer. These Lombok Strait data are adapted from Visser (2004). Computations were conducted using Δx = 300 m, and Δt = 1 s, and absorbing boundaries were applied for both sides. In Fig. 9 snapshots of surfaces (blue lines) and interfaces (red lines) at subsequent times are given. Initially, still water level and interface are located at z = 0 m and z = -150 m, respectively. After the first 3000 s, the interface starts to deform, and it becomes larger as time progresses. These snap shoots describe how internal waves developed as a response of the interaction between horizontal current and the bottom sill. During the whole process, the free surface wave is almost undisturbed. Thus, we have successfully implemented the two-layer non-hydrostatic scheme to simulate the formation of an interfacial wave in response to a current passing through the sill.

Conclusions
The non-hydrostatic model of the two-layer fluid system and its corresponding discrete model have been formulated. Dispersion relation of this scheme has shown to confirm the analytical dispersion relation, over a wide range of wave number. This is validated with a simulation of a baro-   Fig. 9. Snapshots of surface η 1 (x, t) (blue lines) and interface η 2 (x, t) (red lines) at subsequent times. clinic mode oscillation that has a correct period. Then, various nonlinear simulations were performed: the propagation of the internal solitary wave, as well as interface disturbance that interact with a triangular barrier and developed into an interfacial wave.
By Lombok Strait topography with a sill, the interaction of upcoming flow with the sill simulates the formation of an internal wave in Lombok Strait.