A NUMERICAL CALCULATION FOR INTERNAL WAVES OVER TOPOGRAPHY

The internal waves propagating from the deep to shallow, and the shallow to deep, areas in the two-layer fluid systems, have been numerically simulated by solving the set of nonlinear equations, based on the variational principle in consideration of both the strong nonlinearity and strong dispersion of internal waves. The incident wave in the deep area, is the BO-type downward convex internal wave, which is the numerical solution obtained for the present fundamental equations. In the cases where the interface elevation is below, or equal to, the critical level in the shallow area, the disintegration of the internal waves occurs remarkably, leading to a long wave train. The lowest elevation of the interface, increases after its gradual decrease in the shallow area, where the interface is above the critical level, while the lowest elevation of the interface, increases through the internal-wave propagation in the shallow area, where the interface elevation is below, or equal to, the critical level, after its steep decrease around the boundary between the area over the upslope, and the shallow region.


INTRODUCTION
Internal waves show interface deformation when they reach a continental shelf from the deep ocean, as studied by e.g.Ostrovsky & Stepanyants (1989).Grimshaw et al. (2004) numerically simulated the transformation of internal solitary waves using the EKdV eq., obtained by Lee & Beardsley (1974) by extending the KdV eq.However, the applicability of the nonlinear wave equations, derived using perturbation methods for two-layer systems, is restricted owing to the perturbation assumptions, such that, for example, the BO equation (Benjamin, 1967;Ono, 1975), as well as the deeper version from Choi & Camassa (1999), is available for internal waves in offshore deep regions, while the KdV equation, and the shallower version of Choi & Camassa (1999), are adopted for those in shallow water over a continental shelf.Conversely, the equations based on a variational principle, without any assumption concerning wave nonlinearity or dispersion, are expected to be applicable to nonlinear waves propagating from deep to shallow water.In the present study, the set of nonlinear equations based on the variational principle (Kakinuma, 2001), is solved numerically, to simulate internal waves in two-layer systems over topography, in consideration of both the strong nonlinearity, and the strong dispersion, of waves.

FUNDAMENTAL EQUATIONS
The illustration in Fig. 1 is our schematic for the multi-layer system of fluids, represented as i (i = 1, 2, •••, I) from top to bottom, respectively, showing only inviscid and incompressible motion.The thickness of the i-layer in still water, is denoted by hi(x).None of the fluids mix even in motion, and the density in the i-layer, ρi, is spatially uniform and temporally constant in each layer.The fluid motion is assumed to be irrotational, resulting in the existence of velocity potential ϕi, which is expanded into a power series on vertical position z with weightings, in a manner similar to that for the derivation of nonlinear surface-wave equations by Isobe (1995), as (1) where Ni is the number of terms for the expanded velocity potential in the i-layer.
In the i-layer, if both the displacement of one interface, z = ηi,1j (x,t) ( j = 0 or 1), and the pressure on the other interface, pi,j (x,t), are known, then the unknown variables are the velocity potential ϕi (x,z,t) and the interface displacement ηi,j (x,t), where ηi,0 (x,t) and ηi,1 (x,t) are the displacements of the lower and upper interfaces of the i-layer, respectively; pi,0 (x,t) and pi,1 (x,t) are the pressures at the lower and upper interfaces of the i-layer, respectively.Therefore, the definition of the functional for the variational problem in the i-layer, Fi, is as follows (Kakinuma, 2001), by modifying the functional referred to in Luke (1965) for surface waves with the rotational motion of a fluid: (2) (3) where   = ∑ (  −   ) −1 =1 ghk ;  is a partial differential operator in the horizontal x-y plane, i.e., . The gravitational acceleration g is 9.8 m/s 2 .The plane A, which is the orthogonal projection of the object domain on to the x-y plane, is assumed to be independent of time.
If we consider a two-layer fluid system, over a seabed of a fixed profile, i.e., z = η2,0 = b (x), then the Euler-Lagrange equations for the variational problem are reduced to the nonlinear surface/internal-wave equations for two-layer fluid motion as where the surface, and the interface, profiles are described by z = η1,1 = ζ (x,t), and z = η1,0 = η2,1 = η (x,t), respectively; the pressure at the interface is denoted by p = p1,0 = p2,1; z = b (x) is the seabed profile.The sum rule of product is adopted for subscripts αi, i, and i, such that, for example, 1 in the first term on the left-hand side of Eq. ( 4) denotes the power of .After eliminating p from Eqs. ( 6) and ( 8), we obtain (9) In this paper, the top face of the upper layer is assumed to touch a fixed horizontal plate at z = 0, such that  = 0 without any surface wave.In order to represent both the nonlinearity and dispersion of internal waves propagating from deep to shallow, or from shallow to deep, water, the numbers of terms for the expanded velocity potentials introduced in Eq. ( 1), are N1 = 3 and N2 = 5 for the upper and lower layers, respectively, referring to the accuracy of the numerical results obtained by Yamashita & Kakinuma (2015), using the same equations for the internal solitary waves.

NUMERICAL METHOD AND CONDITIONS
The fundamental equations are solved using the finite-difference implicit scheme developed by Nakayama & Kakinuma (2010), to simulate the propagation of internal waves.In the following discussion, the fluid-density ratio ρ1/ρ2 is 0.98. Figure 2 shows the target domain for computation, where the total water depth is uniform in the areas I, III, and V, while the slope gradient is 0.1 and 0.1 in the areas II and IV, respectively.In the deep-water areas I and V, the ratio between the upper-layer thickness h1, and the total water depth h = h1 + h2, i.e., h1/h, is 0.01 in still water.The values of the thickness ratio h1/h in the shallow area III, are shown in Tab. 1.The seabed profile depicted in Fig. 2 is that for Case A, where h1/h = 0.20 in the shallow area III.
The broken line shown in Fig. 2, indicates the elevation of critical level zc, defined by Funakoshi & Oikawa (1986) as (10) In the deep areas I and V for all the present cases, the interface exists above the critical level, such that downward convex internal waves are stable in the areas I and V, while in the shallow area III for Case C, the interface exists below the critical level, such that upward convex internal waves are stable in the area III for Case C. Shown in Fig. 3 is the incident solitary-wave solution, obtained for the same equations, using the numerical method presented by Yamashita and Kakinuma (2015), where the ratio between the upperlayer thickness in still water, and the total water depth, h1/h, is 0.01, and that between the wave height, and the upper-layer thickness in still water, a/h1, is 0.1, in comparison with the corresponding BO solution.The internal-solitary-wave solution, including the weightings fi,αi (x,0) for the expanded velocity potentials, in addition to the interface displacement η (x,0), is given in the deep area I for all the cases.
The lateral boundary condition at x/hI = 0.0 and 80.0, is the perfect-reflection condition.The grid width, and time step, for computation are Δx = 0.01hI, and Δt = 0.01Δx/Ci,0,I, respectively, where Ci,0,I is the phase velocity of linear internal waves in the area I, on the assumption that the internal waves propagate in shallow water.The phase velocity of linear internal waves in shallow water, Ci,0, is evaluated by ( 11)

INTERNAL WAVES PROPAGATING OVER A MOUND
Figure 4 shows the interface profile at each time for Case A, where h1/h = 0.20 in the shallow area III.The internal wave shows disintegration in the shallow area III after passing the slope, resulting in two downward convex internal waves.As the waves propagate over the downslope in the area IV, the wave height of the first internal wave decreases, while that of the second internal wave increases, after which several internal waves, with both a crest and a trough, are generated through disintegration in the deep area V.
Shown in Fig. 5 is the time variation of the interface profile for Case B, where h1/h = 0.50, i.e., h1 = h2, in the shallow area III.The internal wave shows instability in the shallow area III, leading to a wave train, and then the total length of the wave train increases in the areas IV and V.The partial reflection of the internal wave occurs at the boundary between the areas II and III, as well as III and IV, owing to the sudden change in slope gradient.
Figure 6 depicts the time variation of the interface profile for Case C, where h1/h = 0.67 in the shallow area III.The interface shows remarkable disintegration in the shallow area III, following the increase in the wave height of the first internal wave.The partial reflection of the internal wave occurs especially at the boundary between the areas II and III.After the interface elevation becomes below the critical level, the upward convex internal waves are generated in the shallow area III, as shown in the numerical results obtained by Helfrich et al. (1984), as well as the field data from the observation by Orr & Mignerey (2003).
Figure 7 shows the distributions of both the highest elevation ηmax, and the lowest elevation ηmin, of the interfaces for Cases A, B, and C, where h1/h = 0.20, 0.50, and 0.67, respectively, in the shallow area III.In Case A, where the interface exists above the critical level also in the shallow area III, the highest elevation ηmax of the interface remains constant before the internal wave reaches the area IV, which means that the convex direction has been kept downward in the areas from I to III.In Case B, where the stillwater thickness of the upper layer is the same as that of the lower layer, the highest elevation ηmax of the interface increases as the internal waves propagate in the shallow area III, for the wave disintegration  continues owing to the instability in the shallow area III.In Case C, where the interface exists below the critical level in the shallow area III, the remarkable disintegration occurs just after the internal wave passes the boundary between the areas II and III, after which the highest elevation ηmax of the interface decreases as the internal waves propagate in the shallow area III.Conversely, in Case A, the lowest elevation ηmin of the interface increases after its gradual decrease at 14.5 ≤ x/hI ≤ 25.6 in the shallow area III, while in Cases B and C, the lowest elevation ηmin of the interface increases through the internal-wave propagation in the shallow area III, after its steep decrease around the boundary between the areas II and III.

CONCLUSIONS
The internal waves propagating from the deep to shallow, and the shallow to deep, areas in the twolayer fluid systems, were numerically simulated by solving the set of nonlinear equations, based on the variational principle in consideration of both the strong nonlinearity and the strong dispersion of internal waves.The incident wave in the deep area, was the BO-type downward convex internal wave, which was the numerical solution obtained for the present fundamental equations.In the cases where the interface elevation was below, or equal to, the critical level in the shallow area, the disintegration of the internal waves occurred remarkably, leading to a long wave train.The lowest elevation of the interface, increased after its gradual decrease in the shallow area, where the interface was above the critical level, while the lowest elevation of the interface, increased through the internal-wave propagation in the shallow area, where the interface elevation was below, or equal to, the critical level, after its steep decrease around the boundary between the area over the upslope, and the shallow region.

Figure 1 .
Figure 1.A schematic for a multi-layer fluid system.

Figure 2 .
Figure 2. The target domain with a mound in Case A, where hI is the total water depth in the deep areas I and V; h1/h = 0.01 in the deep areas I and V, while h1/h = 0.20 in the shallow area III.The broken line indicates the critical level zc, defined by Eq. (10).

Figure 3 .
Figure 3.The interface profile of the incident internal solitary wave, given in the deep area I shown in Fig. 2, in comparison with the corresponding BO solution, where hI is the total water depth in the deep area I.The fluiddensity ratio ρ1/ρ2 is 0.98; the ratio between the upper-layer thickness in still water, and the total water depth, h1/h, is 0.01; the ratio between the wave height, and the upper-layer thickness in still water, a/h1, is 0.1.

Figure 4 .
Figure 4.The time variation of the interface profile in Case A, where h1/h = 0.20 in the shallow area III shown in Fig. 2; hI is the total water depth in the deep areas I and V; Ci,0,I is the phase velocity of linear internal waves in the deep areas I and V.

Figure 5 .
Figure 5.The time variation of the interface profile in Case B, where h1/h = 0.50 in the shallow area III shown in Fig. 2; hI is the total water depth in the deep areas I and V; Ci,0,I is the phase velocity of linear internal waves in the deep areas I and V.

Figure 6 .
Figure 6.The time variation of the interface profile in Case C, where h1/h = 0.67 in the shallow area III shown in Fig. 2; hI is the total water depth in the deep areas I and V; Ci,0,I is the phase velocity of linear internal waves in the deep areas I and V.

Figure 7 .
Figure 7.The distributions of both the highest elevation ηmax, and the lowest elevation ηmin, of the interfaces in Cases A, B, and C, where h1/h = 0.20, 0.50, and 0.67, respectively, in the shallow area III shown in Fig. 2; hI is the total water depth in the deep areas I and V; the solid and broken lines indicate ηmax and ηmin, respectively.